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The Lorentz covariant theory of propagation of light in the (weak) gravitational fields of N-body 
systems consisting of arbitrarily moving point-like bodies with constant masses m a (a = 1, 2, ..., N) 
is constructed. The theory is based on the Lienard-Wiechert representation of the metric tensor 
which describes a retarded type solution of the gravitational field equations. A new approach for 
integrating the equations of motion of light particles (photons) depending on the retarded time 
argument is invented. Its application in the first post-Minkowskian approximation, which is linear 
with respect to the universal gravitational constant, G, makes it evident that the equations of light 
propagation admit to be integrated straightforwardly by quadratures. Explicit expressions for the 
trajectory of light ray and its tangent vector are obtained in algebraically closed form in terms of 
functionals of the retarded time. General expressions for the relativistic time delay, the angle of 
light deflection, and the gravitational shift of electromagnetic frequency are derived in the form of 
instantaneous functions of the retarded time. They generalize previously known results for the case 
of static or uniformly moving bodies. The most important applications of the theory to relativistic 
astrophysics and astrometry are given. They include a discussion of the velocity dependent terms 
in the gravitational lens equation, the Shapiro time delay in binary pulsars, gravitational Doppler 
shift, and a precise theoretical formulation of the general relativistic algorithms of data processing 
of radio and optical astrometric measurements made in the non- stationary gravitational field of 
the solar system. Finally, proposals for future theoretical work being important for astrophysical 
applications arc formulated. 



I. INTRODUCTION AND SUMMARY 

The exact solution of the problem of propagation of electromagnetic waves in non-stationary gravitational fields 
is extremely important for modern relativistic astrophysics and fundamental astrometry. Until now it are electro- 
magnetic signals coming from various astronomical objects which deliver the most exhaustive and accurate physical 
information about numerous intriguing phenomena going on in the surrounding universe. Present day technology 
has achieved a level at which the extremely high precision of current ground-based radio interferometric astronomical 
observations approaches 1 /iarcsec. This requires a better theoretical treatment of secondary effects in the propagation 
of electromagnetic signals in variable gravitational fields of oscillating and precessing stars, stationary and coalescing 
binary systems, and colliding galaxies fil. Future space astrometric missions like GAIA or SIM J3| will also have 
precision of about 1-10 /uarcsec on positions and parallaxes of stars, and about 1-10 /xarcsec per year for their proper 
motion. At this level of accuracy we are not allowed anymore to treat the gravitational field of the solar system as 
static and spherically symmetric. Rotation and oblateness of the Sun and large planets as well as time variability of 
the gravitational field should be seriously taken into account J4|] . 

As far as we know, all approaches developed for integrating equations of propagation of electromagnetic signals in 
gravitational fields were based on the usage of the post-Newtonian presentation of the metric tensor of the gravitational 
field. It is well-known (see, for instance, the textbooks |§], |§) that the post-Newtonian approximation for the metric 
tensor is valid only within the so-called "near zone" . Hence, the post-Newtonian metric can be used for the calculation 
of light propagation only from the sources lying inside the near zone of a gravitating system of bodies. The near zone 
is restricted by the distance comparable to the wavelength of the gravitational radiation emitted from the system. 
For example, Jupiter orbiting the Sun emits gravitational waves with wavelength of about 0.3 parsecs, and the binary 
pulsar PSR B1913+16 radiates gravitational waves with wavelength of around 4.4 astronomical units. It is obvious 
that the majority of stars, quasars, and other sources of electromagnetic radiation are usually far beyond the boundary 
of the near zone of the gravitating system of bodies and another method of solving the problem of propagation of light 
from these sources to the observer at the Earth should be applied. Unfortunately, such an advanced technique has 
not yet been developed and researches relied upon the post-Newtonian approximation of the metric tensor assuming 
implicitly that perturbations from the gravitational-wave part of the metric are small and may be neglected in the 
data processing algorithms H, @ - H]- However, neither this assumption was ever scrutinized nor the magnitude of 
the neglected residual terms was estimated. An attempt to clarify this question has been undertaken in the paper 
pj where the matching of asymptotics of the internal near zone and external Schwarzschild solutions of equations of 
light propagation in the gravitational field of the solar system has been employed. Nevertheless, a rigorous solution 
of the equations of light propagation being simultaneously valid both far outside and inside the solar system was not 
found. 

One additional problem to be enlightened relates to how to treat the motion of gravitating bodies during the time of 



propagation of light from the point of emission to the point of observation. The post-Newtonian metric of a gravitating 
system of bodies is not static and the bodies move while light is propagating. Usually, it was presupposed that the 
biggest influence on the light ray the body exerts when the photon passes nearest to it. For this reason, coordinates 
of gravitating bodies in the post-Newtonian metric were assumed to be fixed at a specific instant of time t a (see, 
for instance, M - pi]) being close to that of the closest approach of the photon to the body. Nonetheless, it was 
never fully clear how to specify the moment t a precisely and what magnitude of error in the calculation of relativistic 
time delay and/or the light deflection angle one makes if one choses a slightly different moment of time. Previous 
researches gave us different conceivable prescriptions for choosing t a which might be used in practice. Perhaps, the 
most fruitful suggestion was that given by Hellings |1(|] and discussed later on in a paper Q in more detail. This was 
just to accept t a as to be exactly the time of the closest approach of the photon to the gravitating body deflecting 
the light. Klioner & Kopeikin [|4| have shown that such a choice minimizes residual terms in the solution of equation 
of propagation of light rays obtained by the asymptotic matching technique. We note, however, that neither Hellings 
JlO[ nor Klioner & Kopeikin jj] have justified that the choice for t a they made is unique. 

Quite recently we have started the reconsideration of the problem of propagation of light rays in variable gravi- 
tational fields of gravitating system of bodies. First of all, a profound, systematic approach to integration of light 
geodesic equations in arbitrary time-dependent gravitational fields possessing a multipole decomposition [fij , mU has 
been worked out. A special technique of integration of the equation of light propagation with retarded time argument 
has been developed which allowed to discover a rigorous solution of the equations everywhere outside a localized 
source emitting gravitational waves. The present paper continues the elaboration of the technique and makes it clear 
how to construct a Lorentz covariant solution of equations of propagation of light rays both outside and inside of a 
gravitating system of massive point-like particles moving along arbitrary world lines. In finding the solution we used 
the Lienard- Wiechert presentation for the metric tensor which accounts for all possible effects in the description of 
the gravitational field and is valid everywhere outside the world lines of the bodies. The solution, we have found, 
allows to give an unambiguous theoretical prescription for chosing the time t a . In addition, by a straightforward 
calculation we obtain the complete expressions for the angle of light deflection, relativistic time delay (Shapiro effect), 
and gravitational shift of observed electromagnetic frequency of the emitted photons. These expressions are exact at 
the linear approximation with respect to the universal gravitational constant, G, and at arbitrary order of magnitude 
with respect to the parameter v a /c where v a is a characteristic velocity of the a-th light-deflecting body, and c is 
the speed of light [Q. We devote a large part of the paper to the discussion of practical applications of the new 
solution of the equations of light propagation including moving gravitational lenses, timing of binary pulsars, the 
consensus model of very long baseline interferometry, and the relativistic reduction of astrometric observations in the 
solar system. 

The formalism of the present paper can be also used in astrometric experiments for testing alternative scalar- 



tensor theories of gravity after formal replacing in all subsequent formulas the universal gravitational constant G 
by the product G(j* + l)/2, where 7* is the effective light-deflection parameter which is slightly different from its 
weak- field limiting value 7 of the standard parameterized post-Newtonian (PPN) formalism Ej ]. This statement 
is a direct consequence of a conformal invariance of equations of light rays Ji6| and can be immediately proved by 
straightforward calculations. Solar system experiments have not been sensitive enough to detect the difference between 



the two parameters. However, it may play a role in the binary pulsars analysis 14 1. 

The paper is organized as follows. Section 2 presents a short description of the energy-momentum tensor of the 
light-deflecting bodies and the metric tensor given in the form of the Lienard- Wiechert potential. Section 3 is devoted 
to the development of a mathematical technique for integrating equations of propagation of electromagnetic waves in 
the geometric optics approximation. Solution of these equations and relativistic perturbations of a photon trajectory 
are given in Section 4. We briefly outline equations of motion for slowly moving observers and sources of light in 
Section 5. Section 6 deals with a general treatment of observable relativistic effects - the integrated time delay, the 
deflection angle, and gravitational shift of frequency. Particular cases are presented in Section 7. They include the 
Shapiro time delay in binary pulsars, moving gravitational lenses, and general relativistic astrometry in the solar 
system. 

II. ENERGY-MOMENTUM AND METRIC TENSORS 

The tensor of energy-momentum of a system of massive particles is given in covariant form, for example, by Landau 
& Lifshitz 

N 



T^(t,x) = ^f^(i)«5(x-x a (i)) , (1) 

a=l 

ff{t)=m al - l {t)u a a {t)ui{t), (2) 



where t is coordinate time, x = x 1 = (x 1 ,x 2 ,x 3 ) denotes spatial coordinates of a current point in space, m a is the 
constant (relativistic) rest mass of the a-th particle, x a (t) are spatial coordinates of the a-th massive particle which 
depend on time i, v a (t) = dx a (t)/dt is velocity of the a-th particle, "f a (t) = [1 — ^aM] ~^ 2 i s the (time-dependent) 
Lorentz factor, u"(t) — {"f a {t), 7 Q (£)v a (i)} is the four-velocity of the a-th particle, <5(x) is the usual 3-dimensional 
Dirac delta- function. In particular, we have 

rpO0 M m a rf,oi,.s m a v l a (t) -■ m a v l a (t)vj(t) 

\A-W) y/l-vi(t) y/l-VJ(t) 

The metric tensor in the linear approximation reads 



g a p(t, x) = r) a p + h a p(t, x) , (4) 



where r\ a p — diag(— 1, +1, +1, +1) is the Minkowski metric of flat space-time and the metric perturbation h a p(t, x) is 
a function of time and spatial coordinates p8| . It can be found by solving the Einstein field equations which read in 
the first post-Minkowskian approximation and in the harmonic gauge |y| as follows ( |2fj] , chapter 10) 

□ h a p(t, x) = -16irS a p(t, x) , (5) 

where 

S a p(t,x) = T a/3 (i,x) - -rja/3 T\(t,x) . (6) 

The solution of this equations has the form of the Lienard- Wiechert potential [ElJ . In order to see how it looks like 
we represent the tensor of energy-momentum in a form where all time dependence is included in a one-dimensional 
delta-function 

/+oo 
dt'S(t' -t) T Q/3 (t',x) . (7) 

-oo 

Here t' is an independent parameter along the world lines of the particles which does not depend on time t. The 
solution of equation rtq) can be found using the retarded Green function p2[ , and after integration with respect to 
spatial coordinates, using the 3-dimcnsional delta-function, it is given in the form of an one-dimensional retarded-time 
integral 

N ,,+oc 

h^(t,x) = Y. hf(t',t,x)dt' , (8) 

1 J— oo 



ft«"(f,i,x) = 4 



T: p {t')-\^f^{t') 



S[t'-t + r a (t')} 

ra{t>) ' (9j 



where r a (t') = x x Q (t'), and r a (t') — |r a (^')l i s the usual Euclidean length of the vector. 

The integral (pf) can be performed explicitly as described in, e.g., ( |GM, section 14). The result is the retarded 
Lienard- Wiechert tensor potential 

h^(t, x} = 4j2 a ( , 2 tt £b tt > ( 10 ) 

f-[ r a (s) - v (s) • r a {s) 
where the retarded time s — s(t,x) for the a-th body is a solution of the light-cone equation |26| 

s+|x-x a (s)|=t. (11) 



Here it is assumed that the field is measured at time t and at the point x. We shall use this form of the metric 
perturbation h a p(t, x) for the integration of the equations of light geodesies in the next section. It is worth emphasizing 
that the expression for the metric tensor (n(f) is Lorentz-covariant and is valid in any harmonic coordinate system 
admitting a smooth transition to the asymptotically flat space-time at infinity and relating to each other by the Lorcntz 
transformations of theory of special relativity km , p7| . A treatment of post-linear corrections to the Lienard- Wiechert 
potentials ( |l0| ) is given, for example, in a series of papers by Kip Thorne and collaborators p8[ , Q - p3[ . 

III. MATHEMATICAL TECHNIQUE FOR INTEGRATING EQUATIONS OF PROPAGATION OF 

PHOTONS 

We consider the motion of a light particle (photon) in the background gravitational field described by the metric 
(||). No back action of the photon on the gravitational field is assumed. Hence, we are allowed to use equations of 
light geodesies directly applying the metric tensor in question. Let the motion of the photon be defined by fixing the 
mixed initial-boundary conditions (see Fig pi) 

x (to = xo , V - ' = k , (12) 

at 

where k 2 = 1 and, henceforth, the spatial components of vectors are denoted by bold letters. These conditions define 
the coordinates xo of the photon at the moment of emission of light, to, and its velocity at the infinite past and infinite 
distance from the origin of the spatial coordinates (that is, at the, so-called, past null infinity). 

The original equations of propagation of light rays are rather complicated jl| . They can be simplified and reduced 
to the form which will be shown later in this section. In order to integrate them we shall have to resort to a special 
approximation method. In the Minkowskian approximation of the flat space-time the unperturbed trajectory of the 
light ray is a straight line 

x\t) = x%(t) = x l + k* (t - t ) , (13) 



where to, x , and k l = k have been defined in equation (|12[). In this approximation, the coordinate speed of the 
photon is x l = k l and is considered as a constant in the expression for the light-ray-perturbing force. 

It is convenient to introduce a new independent parameter r along the photon's trajectory according to the rule 

r = k-xjv(t) =t-t +k-x , (14) 

where here and in the following the dot symbol " • " between two spatial vectors denotes the Euclidean dot product. 
The time to of the light signal's emission corresponds to the numerical value of the parameter To = k • xo, and the 
numerical value of the parameter t = corresponds to the time 



t* = t - k • x , (15) 

which is the time of the closest approach of the unperturbed trajectory of the photon to the origin of an asymptotically 
flat harmonic coordinate system. We emphasize that the numerical value of the moment t* is constant for a chosen 
trajectory of light ray and depends only on the space-time coordinates of the point of emission of the photon and the 
point of its observation. Thus, we find the relationships 

T = t-t*, Tq = t -t* , (16) 

which reveals that the variable r is negative from the point of emission up to the point of the closest approach 
x l (t*) — £*, and is positive otherwise [|34| . The differential identity dt = dr is valid and, for this reason, the 
integration along the light ray's path with respect to time t can be always replaced by the integration with respect to 
variable r. 

Making use of the parameter t, the equation of the unperturbed trajectory of the light ray can be represented as 

x i (r)=x i N (r) = k i T + i i , (17) 

and the distance, r(r) — |xjv(t)|, of the photon from the origin of the coordinate system reads 

r{r) = \/t 2 + d? . (18) 

The constant vector f = ( = k x (xo x k) = k x (xjv(t) x k) is called the impact parameter of the unperturbed 
trajectrory of the light ray, d — |£| is the length of the impact parameter, and the symbol " x " between two vectors 
denotes the usual Euclidean cross product of two vectors. We note that the vector £ is transverse to the vector k. It 
is worth emphasizing once again that the vector £ l is directed from the origin of the coordinate system towards the 
point of the closest approach of the unperturbed path of the light ray to the origin. This vector plays an auxiliary 
role in our discussion and, in general, has no essential physical meaning as it can be easily changed by the shift of the 
origin of the coordinates |35) . 

Implementing the two new parameters r, £ and introducing the four-dimensional isotropic vector k a = (1, k 1 ) one 
can write the equations of light geodesies as follows (for more details see the paper jl]] and reference | J36| ) 



X\T) = \k a kpd t h a P{T^) - d T 



k a h m (r,0 + ifcWV.C) - ifc J fc P fc^(r,0 



(19) 



where dots over the coordinates denote differentiation with respect to time, d T = d/dr, di = Pijd/d^ , and Pij = 
5%j — kikj is the operator of projection onto the plane being orthogonal to the vector k l , and all quantities on the right 
hand side of equation ([L9) are taken along the light trajectory at the point corresponding to a numerical value of the 
running parameter r while the parameter £ is assumed as constant. Hence, the equation (|19|) should be considered as 



an ordinary, second order differential equation in variable r |37| . The given form of equation ( |19| ) already shows that 
only the first term on the right hand side of it can contribute to the deflection of light if the observer and the source 
of light are at spatial infinity. Indeed, a first integration of the right hand side of the equation dlfl) with respect to 
time from — oo to +00 brings all terms showing time derivatives to zero due to the asymptotic flatness of the metric 
tensor which proves our statement (for more details see the next section). 

However, if the observer and the source of light are located at finite distances from the origin of coordinate system, 
we need to know how to perform the integrals from the metric perturbations (H) with respect to the parameter r 
along the unperturbed trajectory of light ray. Let us denote those integrals as 

# Q/ V,0= f h^[j,iL(p)]da , (20) 

J —OO 



D^( T ,i)= B^(a,£)da, (21) 



where the metric perturbation h a ^[a, x(cr)] is defined by the Lienard- Wiechert potential (|8|) and a is a parameter 



along the light ray having the same meaning as the parameter r in equation (14). In order to calculate the integrals 
(|20|), (|2l|) it is useful to change in the integrands the time argument, er, to the new one, £, defined by the light-cone 
equation (Ww which after substitution for x the unperturbed light trajectory (|17|) reads as follows Eq] 



a + t* = C+IS + kcr-x a (C)|. ( 22 ) 

The differentiation of this equation yields a relationship between differentials of the time variables a and C, an d 
parameters t * , £% k l 

dC, (r a - v a • r a ) = da (r a - k • r a ) + r a dt* - r a • d£ - av a ■ dk , (23) 

where the coordinates, x a , and the velocity, v a , of the a-th body are taken at the retarded time £, and coordinates 
of the photon, x, are taken at the time &(()■ From equation (|23| ) we immediately obtain the partial derivatives with 
respect to the parameters 

cjr 1 

, (24) 

r a - v„ • r a 

and have the relationship between the time differentials along the world line of the photon which reads as follows 

da = d( r a - v ■ r a 
r a - k • r a 

If the parameter a runs from —00 to +00, the new parameter £ runs from C_oo = —00 to (+00 = t* + k • x a (^ +00 ) 
provided the motion of each body is restricted inside a bounded domain of space, like in the case of a binary system. 
In case the bodies move along straight lines with constant velocities, the parameter a runs from —00 to +00, and the 



d( 


r a 




p.. r ] 


d( 


dt* r a 


- v Q • r a ' 


r a - v a • r Q 


dk l 



parameter £ runs from — oo to +00 as well. In addition, we note that when the numerical value of the parameter a 
is equal to the time of observation r, the numerical value of the parameter £ equals to s(t), which is found from the 
equation of the light cone (nil) in which the point x denotes spatial coordinates of observer. 
After transforming time arguments the integrals (EOJ), (Ell) take the form 



B^(s)=J2Bf(s), Bf(s)=A w^JlV^V , (26) 



a=l 



r o (o-,C)-k-r„((T,0 



D aP {s) = Y, BfUa)]da, (27) 



where retarded times in the upper limits of integration depend on the index of each body as it has already been 
mentioned in the previous text. Now we give a remarkable, exact relationship 

r a (a, C) - k • r a (a, C) = t* + k • x a (()-(, (28) 

which can be proved by direct use of the light-cone equation (pf) and the expression (ujh for the unperturbed trajectory 
of light ray. It is important to note that in the given relationship t* is a constant time corresponding to the moment 
of the closest approach of the photon to the origin of coordinate system. The equation (Pq) shows that the integrand 
on the left hand side of the second of equations (£6|) does not depend on the parameter a at all, and the integration 
is performed only with respect to the retarded time variable £. Thus, just as the law of motion of the bodies x (t) 
is known, the integral ( [26| ) can be calculated cither analytically or numerically without solving the complicated light- 
cone equation (0) to establish the relationship between the ordinary and retarded time arguments. This statement 
is not applicable to the integral (g7|) because transformation to the new variable (g5|) does not eliminate from the 
integrand of this integral the explicit dependence on the argument of time t. Fortunately, as it is evident from the 
structure of equation (Ejt), we do not need to calculate this integral. 

Instead of that, we need to know the first spatial derivative of D a P(s) with respect to £\ In order to find it we note 
that the integrand of B a "(s) does not depend on the variable £\ This dependence manifests itself only indirectly 
through the upper limit s(t,£) of the integral because of the structure of the light-cone equation which assumes at 
the point of observation the following form 

r + t* = s+|| + kr-x (s)| . (29) 

For this reason, a straightforward differentiation of B a @(s) with respect to the retarded time s and the implementation 
of formula (g4|) for the calculation of the derivative ds/d^ 1 at the point of observation yields |46| 

b.b>»(,) - if gwzi?£g£w mM oo) 

1 lj ~ ^ r„(a)-kT (*) r (a)-v (*).r (*)- [M) 

10 



This result elucidates that BiB af3 (s) is not an integral but instanteneous function of time and, that it can be calculated 
directly if the motion of the gravitating bodies is given. While calculating diD al3 (s) we use, first, the formula ( pOJ ) 
and, then, replacement of variables (pq). Proceeding in this way we arrive at the result 

„ i •'- v „- ,■• » K(ct,C) -k-r a (o-,C)] 2 



N - T N 



ftD^( a ) = £ / ^B^[C(a)]da = -4£ / a vs 2 ' , aA ^: 2 n- rj(o-,C)dC (31) 



„_,, ., x, [i* + k-x a (C)-C] 2 J J-oo [t* + k-x Q (C)-C] 5 



where the numerical value of the parameter s in the upper limit of the integral is calculated by solving the light-cone 
equation ( |ll| ) . Going back to the equation ( p8| ) we find that the integrand of the integral (|T]) depends only on the 
retarded time argument £. Hence, again, as it has been proven for B a @(s), the integral ( J3l| ) admits a direct calculation 
as soon as the motion of the gravitating bodies is prescribed J17J . 

IV. RELATIVISTIC PERTURBATIONS OF A PHOTON TRAJECTORY 

Perturbations of the trajectory of the photon are found by straightforward integration of the equations of light 
geodesies (Q) using the expressions (|20| ) , ( pi| ) . Performing the calculations we find 

x i ( T )=k i + E i (T), (32) 

x\t)=x* n (t) + Z\t)-Z\t ), (33) 

where r and tq correspond, respectively, to the moment of observation and emission of the photon. The functions 
S*(r) and S l (r) are given as follows 

S l (r) = ^k a k diB^(r) - k a h ai (r) - ^/i 00 (r) + lk%k q h™(r) , (34) 

S 4 (r) = ~k a k diD a P(T) - k a B m {r) - ^k l B 00 ( T ) + l -k l k p k q B™{r) , (35) 



where the functions /i Q/3 (r), B o/j (t), ^""(t), and a,L» a/5 (r) are defined by the relationships ©, @, Q, and 
( pl| ) respectively. 

The latter equation can be used for the formulation of the boundary value problem for the equation of light geodesies. 
In this case the initial position, xo = x(£o), and final position, x = x(i), of the photon are given instead of the initial 
position Xo of the photon and the direction of light propagation k given at past null infinity. All what we need for 
the formulation of the boundary value problem is the relationship between the unit vector k and the unit vector 

|x-x | 
11 



which defines a geometric direction of the light propagation from observer to the source of light in flat space-time (see 
Fig 0). The formulas (§|) and © yield 

fc i = -ir-/3 i (T,0+/3 i (To,|), (37) 

where relativistic corrections to the vector K l are defined as follows 

a 1 (r, i) = 2 |x-x | — ' (38) 

& l (to , i) = ^-^ |x-x | • (39) 

We emphasize that the vectors /3 l (r, £) = /3 and /^(to,^) = /3 are orthogonal to the vector k and are taken at the 
points of observation and emission of the photon respectively. The relationships obtained in this section are used for 
the discussion of observable relativistic effects in the following section. 

V. EQUATIONS OF MOTION FOR MOVING OBSERVERS AND SOURCES OF LIGHT 

The knowledge of trajectory of motion of photons in the gravitational field formed by a N-body system of arbitrary- 
moving point masses is necessary but not enough for the unambiguous physical interpretation of observational effects. 
It also requires to know how observers and sources of light move in the gravitational field of this system. Let us assume 
that observer and the source of light are point-like massless particles which move along time-like geodesic world lines. 
Then, in the post-Minkowskian approximation equations of motion of the particles, assuming no restriction on their 
velocities except for that v < c (see, however, discussion in p0|), read 

**(*) = Tj^oo.i ~ hoi,t - 2 h oo,ti l - h ik , x x k - (h 0i , k - h Qk ^) x k (40) 

- h m ,ki k x l - f h ik j - -h kj A J x k x J + f -h k j,t - h okJ J ^x^x 1 + 0{—) . 



In the given coordinate system for velocities much smaller than the speed of light, the equation (40) reduces to 



x i (t) = -h 00>i + O(-^)+O(-^). (41) 



Regarding specific physical conditions either the post-Minkowskian equation (BOI) or the post-Newtonian equation ( |4l| ) 
should be integrated with respect to time to give the coordinates of an observer, x(t), and a source of light, xo(io)j 
as a function of time of observation, t, and of time of emission of light, £o> respectively. We do not treat this problem 
in the present paper as its solution has been developed with necessary accuracy by a number of previous authors. In 
particular, the post-Minkowskian approach for solving equations of motion of massive particles is thoroughly treated 
in |43|1 , [J44| , |4q| , and references therein. The post-Newtonian approach is outlined in details, for instance, in M, 
||, [|49| - |5l}] , and references therein. In what follows, we assume the motions of observer, x(i), and source of light, 
x o(^o)j to be known with the required precision. 
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VI. OBSERVABLE RELATIVISTIC EFFECTS 
A. Shapiro Time Delay 

The relativistic time delay in propagation of electromagnetic signals passing through the static, spherically- 
symmetric gravitational field of the Sun was discovered by Irwin Shapiro p2j. We shall give in this paragraph 
the generalization of his idea for the case of the propagation of light through the non- stationary gravitational field 
formed by an ensemble of N arbitrary-moving bodies. The result, which we shall obtain, is valid not only when the 
light ray propagates outside the system of the bodies but also when light goes through the system. In this sense we 
extend our calculations made in a previous paper fij] which treated relativistic effects in propagation of light rays only 
outside the gravitating system having a time-dependent quadrupole moment. 

The total time of propagation of an electromagnetic signal from the point Xo to the point x is derived from equations 
(p3|), (p5|). First, we use the equation (|33| ) to express the difference x — xo through the other terms of the equation. 
Then, we multiply this difference by itself using the properties of the Euclidean dot product. Finally, we find the 
total time of propagation of light, t — to, extracting the square root from the product, and using the expansion with 
respect to the relativistic parameter (Gm a )/(c 2 r a ) which is assumed to be small. It results in 

t-t = |x-x |-k-S(r)H-k-S(ro) , (42) 



t-to = |x-xo|+A(t,t ), (43) 

where |x — Xq| is the usual Euclidean distance between the points of emission, xo, and observation, x, of the photon, 
and A(i, to) is the generalized Shapiro time delay produced by the gravitational field of moving bodies 

N 

._, ,, , S.„.k a B a,3 (Tn\ = 2 



1 1 N 

A(Mo) = -k a k p B a "(T) - -k a k B al3 {T O ) = 2^m a B a {s,s ) . (44) 

™ " a— 1 

In the integral 

/• 8 [l-k-v a (Q] 2 d( 
B a {s,s )= — — - j- -, (45) 

Js Vi-^(C) **+k-x a (C)-C 
the retarded time s is obtained by solving the equation (FTA) for the time of observation of the photon, and Sq is found 
by solving the same equation written down for the time of emission of the photon p3| 

s + |x -x a (s )| =*o • (46) 

The relationships (|4Jj), (^4|) for the time delay have been derived with respect to the coordinate time t. The trans- 
formation from the coordinate time to the proper time T of the observer is made by integrating the infinitesimal 
increment of the proper time along the world line x(i) of the observer p7| 
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r = f U - v 2 (t) - hoo[t,x(t)] - 2h 0i [t,x{t)]v i (t) - h l0 [t,x{t)]v\ty(t)\ dt , (47) 

where i; is the initial epoch of observation, and t is a time of observation. 

The calculation of the integral (|4q ) is performed by means of using a new variable 

y = r+k-x a (C)-C, d|, = -[l-k-v (C)]dC. (48) 

so that the above integral pq) reads 

ni \ f 1 - k ' v a(C) rf ( ln ?/) JA , AQ s 

BaM= -L7^m ^ c - (9) 

Integration by parts results in 

B a(s , So ) = -I_L^ ln[r fl ( S )-k.r tt (,)]+ 1 ; k,V ; (8o) ln[r„(»o) - k • r a (* )] (50) 



f ln(r a - k • r a ) . 

/ — ^JrT t k _ v « _ Va x ( k x Vq )1 ' Va d ^ 

Js„ 1 - Wf ) ' 



'-o (i-«2) 

The first and second terms describe the generalized form of the Shapiro time delay for the case of arbitrary moving 
(weakly) gravitating bodies. The last term in the right hand side of (|o|) depends on the body's acceleration and is 
a relativistic correction comparable, in general case, to the main terms of the Shapiro time delay. This correction is 
identically zero if the bodies move along straight lines with constant velocities. Otherwise, we have to know the law 
of motion of the bodies for its calculation. Neglecting all terms of order v\j (? for the Shapiro time delay we obtain 
the simplified expression 

N 



(51) 



A(t, t ) = -2 V>J In r ° ~ , ' r ° - (k ■ v fl ) ln(r a - k • r a ) 

^— ' L r 0a ~ k • r 0a 

a— 1 v 

+ (k • v a0 ) ln(r 0a - k • r 0a ) + f In [t* + k • x a (C) - C] [k • v„(0] d(] , 

where r Q = x x a (s), r 0a = x - x Q (s ), r a = |r a |, r 0a = |r 0a |, v a = x a (s), v a0 = x Q (s ), and the retarded times s 
and So should be calculated from the light-cone equations (O) and (E6J) respectively. The first term on the right hand 
side of the expression (|5l]) for the Shapiro delay was already known long time ago (see. e.g., - H and references 
therein). Our expressions (|50|), (|5l| ) vastly extends previously known results for they are applicable to the case of 
arbitrary-moving bodies whereas the calculations of all previous authors were severely restricted by the assumption 
that either the gravitating bodies are fixed in space or move uniformly with constant velocities. In addition, there was 
no reasonable theoretical prescription for the definition of the body's positions. The rigorous theoretical derivation 
of the formulas ( pO|) and ( |5l| ) has made a significant progress in clarifying this question and proved for the first time 
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that in calculating the Shapiro delay the positions of the gravitating bodies must be taken at the retarded times 
corresponding to the instants of emission and observation of electromagnetic signal. It is interesting to note that in 
the right hand side of (pi]) the terms being linearly dependent on velocities of bodies can be formally obtained in 
the post-Newtonian approximate analysis as well under the assumption that gravitating bodies move uniformly along 



straight lines 54 1 - |56| . We emphasize once again that this assumption works well enough only if the light travel time 
does not exceed the characteristic Keplerian period of the gravitating system. Previous authors were never able to 
prove that the assumption of uniform motion of bodies can be applied, e.g., for treatment of the Shapiro time delay 
in binary pulsars. We discuss this problem more deeply in the next sections of this paper. 

B. Bending of Light and Deflection Angle 

The coordinate direction to the source of light measured at the point of observation x is defined by the four- vector 
p a = (l,p l ), where p l = — x % , or 

p* = _fc<-H*(r,€), (52) 

and where we have put the minus sign to make the vector p % directed from the observer to the source of light. However, 
the coordinate direction p 1 is not a directly observable quantity. A real observable vector towards the source of light, 
s a = (1, s l ), is defined with respect to the local inertial frame of the observer. In this frame s l = —dX % /dT, where T 
is the observer's proper time and X 1 are spatial coordinates of the local inertial frame. We shall assume for simplicity 
that the observer is at rest [£?J with respect to the (global) harmonic coordinate system (t, x l ). Then the infinitesimal 
transformation from {t, x l ) to (T, X 1 ) is given by the formula 

dT = A° dt + A° dx j , dX 1 = A' dt + K\ dx j , (53) 

where the matrix of transformation A"j depends on the space-time coordinates of the point of observation and is 
defined by the requirement of orthonormality 

g a p = v^K^ ■ (54) 



In particular, the orthonormality condition (54) pre-assumes that spatial angles and lengths at the point of observations 
are measured with the help of the Euclidean metric Sij. For this reason, as the vector s a is isotropic, we conclude 
that the Euclidean length |s| of the vector s 1 is equal to 1. Indeed, one has 

7 ?Q/3 s Q s /3 = -1 + s 2 = 0. (55) 

Hence, |s| = 1, and the vector s points out the astrometric position of the source of light on the unit celestial sphere 
attached to the point of observation. 
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In linear approximation with respect to G the matrix of transformation is as follows pj] 

A° = l- -hoo(t,x) , 
A°. = -/*»(*, x), 

A* = , 



A i i = ^ + -/i ij -(t,x). (56) 



Using the transformation (53) we obtain the relationship between the observable unit vector s l and the coordinate 
direction p % 

A\ pi - A' 

In linear approximation it takes the form 

a* = (l+ h 00 - h 0j A p* + Ihirf . (58) 

Remembering that |s| = 1, we obtain for the Euclidean norm of the vector p l 

IPl = 1 - 2^°o + h Q] p> - -h llP l p> , (59) 

which brings equation (Bq) to the form fl58| 

s i = m l + -P l: >m q h jq (t 7 x) , (60) 

with the Euclidean unit vector m % — p J /|p|- 

Let now denote by a 1 the dimensionless vector describing the angle of total deflection of the light ray measured at 
the point of observation and calculated with respect to vector k l given at past null infinity. It is defined according to 
the relationship 0j 

a l (T ; |) = fc l [k-H(r,0]-S 4 (T,0, (61) 

or 

a i (T,i) = -P i j &(r,0. (62) 

As a consequence of the definitions (J52|) and (|6^) we conclude that 

m l = -k l + a l ( T ,$,) . (63) 

Taking into account expressions (J57|), (p9|), (p2[), and ( p7\ ) we obtain for the observed direction to the source of light 
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s* (r,i)=K i + a 1 (r, i) + P {t,Z)-P\t g ,£)+ f (r, £) , (64) 

where the relativistic corrections /3 l are defined by the equation (p8|) and where 



7 < (r,€) = -^'*^i,(*,x) (65) 

describes the light deflection caused by the deformation of space at the point of observations. If two sources of light 
are observed along the directions s\ and s\, correspondingly, the measured angle tp between them is defined in the 
local inertial frame as follows 

cost/' = si • s 2 , (66) 

where the dot denotes the usual Euclidean scalar product. It is worth emphasizing that the observed direction to 
the source of light (M) includes the relativistic deflection of the light ray which depends not only on quantities taken 
at the point of observation but also on those /3 1 (tq,£) taken at the point of emission of light. Usually this term is 
rather small and can be neglected. However, it becomes important in the problem of propagation of light in the field 
of gravitational waves [Qj or for a proper treatment of high-precision astrometric observations of objects being within 
the boundary of the solar system. 

Without going into further details of the observational procedure we, first of all, give an explicit expression for the 
angle a l (r) 

a\r) = ~k a kpdiB a \T) + k a P) h a i(r) . (67) 

The relationships (jTu), (BO) along with the definition of the tensor of energy-momentum (H) allow to recast the 
previous expression into the form 

a - ( ^ 2 y ,"» ( 1 - k - v -) 2 PjA - 4 y m - 1 ~ k - v ° p\vi, (68) 

where all the quantities describing the motion of the a-th body have to be taken at the retarded time s which relates 



to r = t — t* by the light-cone equation (11). Neglecting all terms of the order v a /c we obtain a simplified form of 
the previous expression 






)= 2 E— 7 VtT> (69) 

^-J r a (r a - k • r ) 

a— 1 v ' 

which may be compared to the analogous expression for the deflection angle obtained previously by many other 
authors in the framework of the post-Newtonian approximation (see [|3| , and references therein) . We note that all 
previous authors fixed the moment of time, at which the coordinates x a of the gravitating bodies were to be calculated 
rather arbitrarily, without having rigorous justification for their choice. Our approach gives a unique answer to this 
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question and makes it obvious that the coordinates x a should be fixed at the moment of retarded time s relating to 
the time of observation t by the light-cone equation (DJ 

The next step in finding the explicit expression for the observed coordinate direction s l is the computation of the 
quantity (3 l (r) given in (pq). We have from formulas d29), (|3l|) the following result for the numerator of (3 1 (t) 



1 



JV 



A' 



-k a k d z D a P(T) - k a P)B aj (T) = -2j2 m " £ iG, »(») - P ) D i( s ) +4^m Q P)El(s) , 



a=l 



where the integrals C a (s), D 3 a (s) and El(s) read as follows 

l-k-v„(C) 



C a {s) 



t* + k-x Q (C)-C 



<K 



0=1 



>/i-«2(0' 



(70) 



(71) 



^iw 



l-k-v (C) 



t* + k-x a (C)-C 



4(0 



Vi - «2(0 



dC 



(72) 



*£00 



l-k-v„(C) 



«i(0 



t* + k-x a (C)-CVi-^(0 



d( 



(73) 



Making use of the new variable y introduced in ( 48 ) and integrating by parts yields 
1 1 - k • v a , f s [k - v a - v a x (k x v a )] • v a d( 



C a (s) 



Di(s) 



y/l - v 2 a r a . - k • r a 



r a - k • r a 



ft 2\3/2 ' 



1 - k • v„ 



y/l - v\ r a - k • r a 



[k - v Q - v Q x (k x v Q )] ■ v a x 3 a d( 



r a -k-r a 



El(s) = —7=2=3 ln(r a - k • r a ) + / ln(r a - k • r„) n£ t, a fc 

V 1 - V a J -co (1 ~ «2) 



V n, k 



(1-vlf 2 



<K 



Ei{s) , 



2\V2 ' 



(74) 



(75) 



(76) 



where Ilfc(C) = ^i- + u *(C) u fc(C) is the spatial part of the operator of projection onto the plane being perpendicular 
to the world line of the a-th body, and the bodies' coordinates and velocities in all terms, being outside the signs of 
integral, are taken at the moment of the retarded time s. The equations (|74|)-(|7q) will be used in section 7 for the 
discussion of the gravitational lens equation with taking into account the velocity of the body deflecting the light rays. 
Finally, the quantity 7 1 (t) can be explicitly given by the following expression 



A 



^ m (k • v ) (P*. vi) 

( T ) = - 2 Z^ R 72 

rr? \/i - vi 



(77) 



where coordinates and velocities of the bodies must be taken at the retarded time s according to equation ( |ll| ) . We 
note that 7* is a very small quantity being proportional to the product (Gm a / 'c 2 r a ){v a / c) . 
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C. Gravitational Shift of Frequency 

The exact calculation of the gravitational shift of electromagnetic frequency between emitted and observed photons 
plays a crucial role for the adequate interpretation of measurements of radial velocities of astronomical objects, 
anisotropy of electromagnetic cosmic background radiation (CMB), and other spectral astronomical investigations. 
In the last several years, for instance, radial velocity measuring technique has reached unprecedented accuracy and is 
approaching to the precision of about 10 cm/sec [ p9| . In the near future there is a hope to improve the accuracy up 
to 1 cm/sec [pOj when measurement of the post-Newtonian relativistic effects in optical binary and/or multiple star 
systems will be possible [ pl| . 

Let a source of light move with respect to the harmonic coordinate system (t, x l ) with velocity Vo(to) = dxo(to)/dto 
and emit electromagnetic radiation with frequency Vq — l/(<57o), where to and Tq are coordinate time and proper 
time of the source of light, respectively. We denote by v — 1/(ST) the observed frequency of the electromagnetic 
radiation measured at the point of observation by an observer moving with velocity v(t) = dx./dt with respect to the 
harmonic coordinate system (£, x l ). We can consider the increments STq and ST as infinitesimally small. Therefore, 
the observed gravitational shift of frequency 1 + z — v jvo can be defined through the consecutive differentiation of 
the proper time of the source of light, Tq, with respect to the proper time of the observer, T, JO - [R4J 

dT dto dt dT 



where the derivative 
dt 



1/2 

1 - v$(t ) - froo(*o,x ) - 2/i 0l (£ ,x ) v l {t ) - hij(t ,xo) v l {t ) v 3 (t Q ) , (79) 



is taken at the point of emission of light, and the derivative 

^= [l-v 2 (t)-h 00 (t,x)-2h 0i (t^)v i (t)-h ij (t,K)v i (t)v^t)y 1/2 , (80) 

is calculated at the point of observation. 

The time derivative along the light-ray trajectory is calculated from the equation ( J43J ) where we have to take into 
account that the function B a (s,So) depends on times t a and t not only through the retarded times Sq and s in the 
upper and lower limits of the integral ( p5[ ) but through the time t* and the vector k both being considered in its 
integrand as time-dependent parameters. Indeed, the infinitesimal increment of times to and/or t causes variations in 
the positions of the source of light and/or observer and, consequently, to the corresponding change in the trajectory 
of light ray, that is in t* and k. Hence, the derivative along the light ray reads as follows 



1 + K • v - 2^ m a 

dto a=l 



x r ds d ds d dt* d dk { d 



dt ds dt dso dt dt* dt dU 



dt N 



1 + K • v + 2^ m a 



ds d dso d dt* d dk l d 
dto ds dto dso dto dt* dto dk l 
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B a {s,so,t*,'k) 

, (81) 



B a (s,so-,t*,\s) 



where the unit vector K is defined in (J3q) and where we explicitly show the dependence of function B a on all 



parameters which implicitly depend on time |65|. 



The time derivative of the vector k is calculated using the approximation k = — K and formula ( pq ) where the 
coordinates of the source of light, xo(to), and of the observer, x(£), are functions of time. It holds 

dk l _ (k x (v x k))' dti _ (k x (vq x k)Y 

~dt~ R ' dh~ R ' ( ' 

where R = |x — xo| is the distance between the observer and the source of light. The derivatives of retarded times s 
and so with respect to t and to are calculated from the formulas ( |ll|) and (|46| ) where we have to take into account 
that the spatial position of the point of observation is connected to the point of emission of light by the unperturbed 
trajectory of light, x(t) = Xo(to) + k (£ — to)- More explicitly, we use for the calculations the following relationships 
@ 

s + |xo(*o)+k(*j*o) (*-*o)-x a (s)| =* , and s + |x (i ) - x o (s )| = *o , (83) 

where the unit vector k must be considered as a two-point function of times t, t with derivatives being taken from 



(p2|). The physical meaning of relationships (p3J) and (277) is the preservation of the intersection at the point of 



observation x(i) of two of the lines forming light cones which relate to propagation of the gravitational field and 
electromagnetic signals, and having vertices at points x n (s) and Xo(to), respectively. Calculation of infinitesimal 
variations of equations (B3[) immediately gives 

ds r Q -k-r a (k x v) • (k x r a ) 
ot r a - v Q • r a r a - v a • r a 

ds_ = (l-k-v )(k-r n ) 
dt r a -v a -r a 

ds a r 0a - v • r 0a -.„. 

-zr = > ( 86 ) 

oto r Qa - v a0 • r 0a 

Time derivatives of the parameter t* are calculated from its original definition t* = to — k • Xo(to), which naturally 
appears in integrands of all integrals, and read 

^ = 1 - k ' V0 + ^T' eT = — r-' (88) 

where the terms of order £/R in both formulas relate to the time derivatives of the vector k . 
Partial derivatives of the function B a (s, sq, t*, k) defined by the integral (|45|) read as follows 
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ds yi^I r a -k-r a ' l8Jj 

dB a 1 (l-k-v a0 ) 2 , , 



d3g 

dB a 

The partial derivative dB a /dt* is found with the help of relationships (|7l|), (|74|). Calculation of the partial derivative 
8B a /dk l is realized by making use of (|72), (73) and (|7q), ( |7q ) respectively. The integrals in (|7l|)-(|76|) are not calculable 



C a (s ) - C a (s) , (91) 



= Di(s )-D i a (s) + 2[Ei( So )-Ei( S )] . (92) 



analytically in general. If we assume that the accelerations of gravitating bodies are small so that the velocity of each 
body can be considered as a constant, the derivatives (|9l|), (B2h are approximated by simpler expressions 

dB a 1 l-k-v„ 1 l-k-v o0 



dt* tJI-vI r a - k • r a ^/i - V 2 Q r 0a - k • r 0a 

dB a 1 - k ■ v a x J a (s) 1 - k • v o0 x J a (s ) 



dk l yjl — v* r a - k • r a ^J\ - V ^ Q r 0a - k • r 0a 



2yi 2v 3 

: ln(r a - k • r a ) a0 ln(r 0a - k • r 0a ) 



(93) 



(94) 



Residual terms, denoted by ellipses, can be calculated from the integrals in (f74[)-(|76|) if one knows the explicit functional 
dependence of the bodies' velocities on time. One expects the magnitude of the residual term to be so small that it 
is unimportant for the following discussion [p7|. The expressions (p3), (M) will be explicitly used in section VII. B for 
discussion of the gravitational shift of frequency by a moving gravitational lens. 

VII. APPLICATIONS TO RELATIVISTIC ASTROPHYSICS AND ASTROMETRY 
A. Shapiro Time Delay in Binary Pulsars 

1. Approximation Scheme for Calculation of the Effect 

Timing of binary pulsars is one of the most important methods of testing General Relativity in the strong grav- 
itational field regime ( [pq| - u% , and references therein) . Such an opportunity exists because of the possibility to 
measure in some binary pulsars the, so-called, post-Keplerian (PK) parameters of the pulsar's orbital motion. The 
PK parameters quantify different relativistic effects and can be analyzed using a theory-independent procedure in 
which the masses of the two stars are the only dynamic unknowns |73J . Each of the PK parameters depends on the 
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masses of orbiting stars in a different functional way. Consequently, if three or more PK parameters can be measured, 
the overdetermined system of the equations can be used to test the gravitational theory. 

Especially important for this test are binary pulsars on relativistic orbits visible nearly edge-on. In such systems 
observers can easily determine masses of orbiting stars measuring the " range" and " shape" of the Shapiro time delay in 
the propagation of the radio pulses from the pulsar to the observer independently of other relativistic effects. Perhaps, 
the most famous examples of the nearly edge-on binary pulsars are PSR B1855+09 and PSR B1534+12. The sine of 
inclination angle, i, of the orbit of PSR B1855+09 to the line of sight makes up a value of about 0.9992 and the range 
parameter of the Shapiro effect reaches 1.27 /is [[74|. The corresponding quantities for PSR B1534+12 are 0.982 and 
6.7 fj,s |75| . All binary pulsars emit gravitational waves, a fact which was confirmed with the precision of about 0.3% 



by Joe Taylor and collaborators |76|. New achievements in technological development and continuous upgrading the 
largest radio telescopes extend our potential to measure with a higher precision the static part of the gravitational 
field of the binary system as well as the influence of the velocity-dependent terms in the metric tensor, generated by 
the moving stars, on propagation of radio signals from the pulsar to the observer. These terms produce an additional 
effect in timing observations which will reveal itself as a small excess to the range and shape of the known Shapiro 
delay making its representation more intricative. The effect under discussion can not be investigated thoroughly 
and self-consistently within the post-Newtonian approximation (PNA) scheme even if the velocity-dependent terms 



in the metric tensor are taken into account |54J] - |56|. This is because the PNA scheme does not treat properly all 
retardation effects in the propagation of the gravitational field. 

In this section we present the exact Lorentz covariant theory of the Shapiro effect which includes, besides of the 
well known logarithm, all corrections for the velocities of the pulsar and its companion. However, later on we shall 
restrict ourselves to terms which are linear with respect to the velocities. The matter is that due to the validity of the 
virial theorem in gravitational bound systems the terms being quadratic with respect to velocities are proportional 
to the gravitational potential of the system. It means that the proper treatment of quadratic with respect to velocity 
terms can be achieved only within the second post-Minkowskian approximation for the metric tensor which is not 
considered in the present paper. 

The original idea of the derivation of the relativistic time delay in the static and spherically symmetric field of a 
self-gravitating body belongs to Irwin Shapiro ]52|] . Regarding binary pulsars the static part of the Shapiro time delay 
has been computed by Blandford & Tcukolsky [J77| under the assumption of everywhere weak and static gravitational 
fields. Nordtvedt |54| ], Klioner pq|, and Wex p& calculated the Shapiro time delay in the gravitational field of 
uniformly moving bodies but without accounting for the retardation in the propagation of the gravitational field. The 
mathematical technique of the present paper allows to treat the relativistic time delay rigorously and account for 
all effects caused by the non- stationary part of the gravitational field of a binary pulsar, that is to find in the first 
post-Minkowskian approximation all special-relativistic corrections of order v a /c, v 2 /c 2 , etc. to the static part of the 
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Shapiro effect where v a denotes characteristic velocity of bodies in the binary pulsar. 

Let us assume that the origin of the coordinate system is at the barycenter of the binary pulsar. Radio pulses are 
emitted rather close to the surface of the pulsar and the coordinates of the point of emission, Xo, can be given by the 
equation 

x =x p (t ) + X(i ), (95) 

where x p are the barycentric coordinates of the pulsar's center-of -mass, and X are the barycentric coordinates of 
the point of emission both taken at the moment of emission of the radio pulse, to- At the moment of emission the 
spatial orientation of the pulsar's radio beam is almost the same with respect to observer at the Earth. Hence, we 
are allowed to assume that the vector X is constant at every "time" when an emission of a radio pulse takes place 
|fr8|| . In what follows the formula ( |5l|) plays the key role. However, before performing the integral in this formula it is 
useful to derive the relationship between retarded times s and so given by the expressions ( |ll| ) and ( [46| ) respectively. 
Subtracting the equation ( jig ) from (^TJ) and taking into account the relationship (f43|), we obtain 

s-s a = R-r a +r Qa + A(t,t a ) , (96) 

where R — |R|, R = x Xo, r a = |x — x a (s)|, and r^ a = |xo — x a (so)|- We note that the point of observation, x, 
is separated from the binary system by a very large distance approximately equal to R. On the other hand, the size 
of the binary system can not exceed the distance ro a . Thus, the Taylor expansion of r a with respect to the small 
parameter ro a /R is admissible. It yields 

r a = |R + x -x a (.s)|=i?-K-[x -x a (.s)]+o(^) , (97) 



where the unit vector K is defined in (|36|). Using the approximation K = k + O(G), formula fl9q ) is reduced to the 
form 

s - s = r 0a - k ■ r 0a + k ■ [x a (s) - x o (s )] + O (^) + 0(G) , (98) 

which explicitly shows that the difference between the retarded times s and Sq is of the order of time interval being 
required for light to cross the binary system. It is this interval which is characteristic in the problem of propagation 
of light rays from the binary (or any other gravitationally bound) system to the observer at the Earth. Therefore, the 
retarded time s taken along the light ray trajectory changes only a little during the entire process of propagation of 
light from the pulsar to the observer while the coordinate time t changes enormously. This remarkable fact was never 
noted in any of previous works devoted to study of propagation of electromagnetic signals from remote astronomical 
systems to observer at the Earth. 

In addition to the expression (g8|) , we can show that time differences sq — to and s — to are also of the same order 
of magnitude as s — So- Indeed, assuming that the velocities of pulsar and its companion are small compared to the 
speed of light, we get from fllq ) and (^8|) for these increments 
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■so - t = -|x - x Q (s )| = -p 0a -p 0a -v a + O (v 2 a ) + 0(G) , (99) 

s-to = -(k ■ p 0o )(l - k • O + K) + 0(G) , (100) 



where p 0a = xo — x a (to), po a = |Poa|i an< ^ ^a — v a(to)- The relationships (fey), (JlOOl) prove our previous statement 
and reveal that coordinates of bodies comprising the system and their time derivatives can be expanded in Taylor 
series around the time of emission of the radio signal to in powers of s — to and/or sq — to. Fig. H illustrates geometry 
of the mutual positions of the binary pulsar and the observer and Fig. explains relationships between position of 
photon on the light trajectory and retarded positions of pulsar and its companion. 

In what follows we concentrate our efforts on the derivation of the linear with respect to velocity of moving bodies 
corrections to the static part of the Shapiro delay. Calculations are realized using the expression (|5l| ) where the 
integral is already proportional to the ratio v a /c. Hence, in order to perform the integration we take into account 
only first terms in the expansion of the integrand with respect to to- Then, the integral reads as 

f \n(r a -k-r a )(k-v a )d( = k-v a (to) I ln[t* + k • x a (t ) - Q d( . (101) 

J S J 8 

After this transformation the integral acquires table form and its calculation is rather trivial. Accounting for (Bq)- 



(100), the result of integration yields 

f In [t* + k • x a (< ) - C] d( = (r 0a - k • r 0a ) [ln(r 0a - k • r 0a ) - 1] - (r„ - k • r a ) ln(r a - k ■ r«) + 0(v a ) , (102) 



where r a and ro a have the same meaning as in (B6|). The result (102) is multiplied by the radial acceleration of 



the gravitating body according to (101). Terms forming such a product can reach in a binary pulsar the maximal 
magnitude of order (Gm a / c 3 )(x / Pb)(v / c) ln(l — sinz), where x is the projected semimajor axis of the binary system 
expressed in light seconds, Pi, is its orbital period, and i is the angle of inclination of the orbital plane of the binary 
system to the line of sight. For a binary pulsar like PSR B1534+12 the terms under discussion are about 10~ 5 /is 
which is too small to be measured. For this reason, all terms depending on the acceleration of the pulsar and its 
companion will be omitted from the following considerations. 

Let us note that coordinates of the a-th body taken at the retarded time s can be expanded in Taylor series in the 
neighborhood of time so 

x Q (s) = x a (s ) + v a (s )(s - so) + 0[(s - s Q ) 2 } , (103) 

or, accounting for (Bq), 

X a ( S ) = X a ( So ) + V a ( P 0a - k • Po a ) + 0(v\) . (104) 

Making use of this expansion one can prove that the large distance, r , relates to the small one, ro a , by the important 
relationship 
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r 2 a - (k • v a f = r 2 a - (k • r 0a ) 2 - 2(p 0a - k ■ p 0a ) [ Va • p 0a - (k ■ p J(k ■ «„)] + 0(t£) . (105) 

Moreover, 

r 0a + k • r 0a = poa + k • p 0a + V a • P 0a + ( k ' ^a)POa + Ofaa) . (106) 

As a concequence of simple algebra we obtain 



kr a r\ -(k-r a ) 2 r 0a + k • r 0a 



r 0Q -k-r 0a r2 Q -(k-r 0a )- r a + k • r Q 



(107) 



which gives after making use of ( |105| ), (106) the following result 
r a - k • r a 1 + k • v a 



r 0a - k • r 0a r a + k ■ r a 
It is straightforward to prove that 



POa + k • p 0a - U a • p 0a + (k • p J(k • Va )] + O (<) . (108) 



r a +k-r a = 2(R + k-r aa ) + Ol -^1 , (109) 

where i? = |R| is the distance from the point of emission to the point of observation. This distance is expanded as 

R = 7^ + X£;+w-Xp-X, (110) 

where TZ is the distance between the barycenters of the binary pulsar and the solar system, x^; is the distance from the 
barycenter of the solar system to the center of mass of the Earth, w is the geocentric position of the radio telescope, 
x p are coordinates of the center of mass of the pulsar with respect to the barycenter of the binary system, and X are 
coordinates of the point of emission of radio pulses with respect to the pulsar proper reference frame. The distance TL 
is gradually changing because of the proper motion of the binary system in the sky. It is well known that the proper 
motion of any star is small and, hence, can be neglected in the time delay relativistic corrections. All other distances 



in formula (110) are of order of either diurnal, or annual, or pulsar's orbital parallax with respect to the distance 7?.. 
Hence, when considering relativistic corrections in the Shapiro time delay, the distance R. can be taken as a constant. 
Such an approximation is more than enough to put 

ln(r„ + k • r„) - ln(2ft) + O (^) ~ const. , (111) 

where 1Z = |72.|. Constant terms are not directly observable in pulsar timing because they are absorbed in the initial 
rotational phase of the pulsar. For this reason, we shall omit for simplicity the term ln(r a + k • r a ) from the final 
expression for the Shapiro time delay. 



Accounting for all approximations having been developed in this section we obtain from (pi]), (108), and (111) 
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A(Mo) = -2^mJ (1 - k ■ v a )ln[p 0a + k ■ p 0a - v a ■ p 0a + (k ■ p 0a ){k ■ v a )} + k ■ v a L 

a=l *■ > 



(112) 



-o ( Grna V *\ + o ( Grria Va x \ + o ( Grria x 

\ c 3 c PbJ \ c 3 1Z 
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This formula completes our analytic derivation of the velocity-dependent corrections to the Shapiro time delay in 
binary systems. It also includes residual terms which have not been deduced by other authors |84[ ]. 

2. Post-Newtonian versus post-Minkowskian calculations of the Shapiro time delay in binary systems 

Our approach clarifies the principal question why the post-Newtonian approximation was efficient for the correct 



calculation of the main (velocity-independent) term in the formula (112) for the Shapiro time delay in binary systems. 
We recall that the post-Newtonian theory operates with the instantaneous values of the gravitational potentials in the 
near zone of the gravitating system. In the post-Newtonian scheme coordinates and velocities of gravitating bodies, 
being arguments of the metric tensor, depend on the coordinate time t. Thus, if we expand these coordinates and 
velocities around the time of emission of light, io, we get for the components of metric tensor g a f3[t, x(i), x (£), v a (i)] 
a Taylor expansion which reads as follows 

g a p[t,x(t),x a (t),v a {t)] = g a /3[t,x(t),x.a(to),v a (to)] + (113) 

dg a p[t,x(t),x a (t ),v a (to)} t dg a0 [t,x(t),x a (to),v a (t o )} if A, , 

Q^ ««(*0) + Q^ ".(to) j(« - *°) + - ■ 

This expansion is divergent if the time interval t — to exceeds the orbital period P of the gravitating system. This is the 
reason why the post-Newtonian scheme does not work if the time of integration of the equations of light propagation 
is bigger than the orbital period. 

On the other hand, the post-Minkowskian scheme gives components of the metric tensor in terms of the Lienard- 
Wiechert potentials being functions of retarded time s. We have shown that in terms of the retarded time argument 
the characteristic time for the process of propagation of light rays from the pulsar to observer corresponds to the 
interval of time being required for light to cross the system. During this time gravitational potentials can not change 
their numerical values too much because of the slow motion of the gravitating bodies. Hence, if we expand coordinates 
of the bodies around to we get for the metric tensor expressed in terms of the Lienard- Wiechert potentials the following 
expansion 

g a p[t, x(t), x a (s), v a (s)] = g a p[t, x(i), x a (i ), v a (i )] + (114) 



dg a/3 [t,x(t),TCa(to),Va(to)] if , dg a p[t, x(f), X a (f ), Vq(to)] .j,. -,(, , , 

v a{ t Q) H KJ Va^O) >(S - to) + 



dxi ay "' dvi 

which always converges because the time difference s — to never exceeds the orbital period (see equation (|100|)). 



Nevertheless, as one can easily see, the leading terms in the expansions ( 113 ) and (114) coincide exactly which 
indicates that the terms in the solution of the equations of light propagation depending only on the static part of 
gravitational field should be identical independently on what kind of approximation scheme is used for finding the 

2G 



metric tensor. Thus, the post-Newtonian approximation works fairly well for finding the leading part of the solution 
of the equations of light geodesies. However, it can not be used for taking into account perturbations of the light 
trajectory caused by the motion of massive bodies in the light-deflecting, gravitationally bounded astronomical systems 

©■ 

It is worth emphasizing once again that our approach is based on the post-Minkowskian approximation scheme for 
the calculation of gravitational potentials which properly accounts for all retardation effects in the motion of bodies 
by means of the Lienard-Wiechert potentials. 

3. Shapiro Effect in the Parametrized Post-Keplerian Formalism 

The parametrized post-Keplerian (PPK) formalism was introduced by Damour & Deruelle [p5| and partially im- 



proved by Damour & Taylor 73 1 . It parametrizes the timing formula for binary pulsars in a general phenomenological 
way [ p7[ . In order to update the PPK presentation of the Shapiro delay we use expression ( p_12[ ) . A binary pulsar 
consists of two bodies - the pulsar (subindex "p") and its companion (subindex "c"). The emission of a radio pulse 
takes place very near to the surface of the pulsar and, according to (pq) and the related discussion, we can approx- 
imate X = Ak where X is the distance from the center of mass of the pulsar to the pulse-emitting point. In this 
approximation we get p Q — Ak and, as a consequence, 

In [p Qp + k ■ p 0p - v p ■ p 0p + (k ■ p 0p )(k ■ v p )] = ln(2A) = const. . (115) 



Hence, the formula (112) for the Shapiro time delay can be displayed in the form 

A(i, t ) = -2mJ (1 - k • v c ) In [p 0c + k • p 0c - v c ■ p 0c + (k • p 0c )(k • v c )] + k • v c I (116) 

-2m p [(1 - k • v p ) ln(2A) + k • v p ] , 

where we have omitted residual terms for simplicity. It was shown in the paper [p9fl that any constant term multiplied 
by the dot product k • v p or k • v c is absorbed into the epoch of the first pulsar's passage through the periastron. 



Thus, we conclude that terms relating to the pulsar in the formula (116) and the very last term in the curl brackets 
are not directly observable. For this reason, we shall omit them in what follows and consider only the logarithmic 
contribution to the Shapiro effect caused by the pulsar's companion. According to formula ( |95| ) we have 

V 

p 0c = r + Ak , POc = r + — k-r + ... , (117) 

r 

where r = x p (to) — x c (to) is the vector of relative position of the pulsar with respect to its companion, r = |r|, and dots 
denote residual terms of higher order. Taking into account all previous remarks and omitting directly unobservable 
terms we conclude that the Shapiro delay assumes the form 



A(i, t ) = -2m c (l - k • v c ) In 



(H I (r + k • r) - v c ■ r + (k • r)(k ■ v c 

r 
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(118) 



If the pulsar's orbit is not nearly edgewise and the ratio X/r is negligibly small the time delay can be decomposed 
into three terms 

A(t, t ) = -2m c In (r + k • r) + 2m c (k ■ v c ) In (r + k • r) + 2m c — V - n ^ . (119) 

r + k ■ r 



The first term on the right hand side of (119) is the standard expression for the Shapiro time delay. The second 
and third terms on the right hand side were discovered by Nordtvedt [B4J and Wex Bq] under the assumption of 
uniform and rectilinear motion of pulsar and companion in the expression for the post-Newtonian metric tensor of the 
binary system. One understands now that this assumption was equivalent to taking into account primary terms of 
retardation effects in propagation of gravitational field of pulsar and its companion. Nevertheless, the approximation 
used by Nordtvedt and Wex works fairly well only for terms linear with respect to velocities of bodies. Had one tried 
to take into account quadratic terms with respect to velocities using the post-Newtonian approach an inconsistent 
result would have been obtained, at least under certain circumstances | p0| . 

In what follows only the case of the elliptic motion of the pulsar with respect to its companion is of importance. 



Moreover, we do not use the expansion (119) keeping in mind the case of the nearly edgewise orbits for which the 
magnitude of r+k-r term can be pretty small near the event of the superior conjunction of pulsar and companion. The 
size and the shape of an elliptic orbit of the pulsar with respect to its companion are characterized by the semi-major 
axis an and the eccentricity e (0 < e < 1). The orientation in space of the plane of the pulsar's motion is defined with 
respect to the plane of the sky by the inclination angle i and the longitude of the ascending node £2. For orientation of 
the pulsar's position in the plane of motion one uses the argument of the pericenter to. More precisely, the orientation 
of the orbit is defined by three unit vectors (1, m, n) having coordinates 0, [p5| 

1 = (cos fl, sinfi, 0) , 
m = (— cos i sin fi, cos i cos $7, sin i) , (120) 

n = (situ sin £7, — sini cos $7, cosi) . 

In this coordinate system we have the unit vector k to be k = K = (0, 0, —1) |9l]]. The coordinates of the pulsar in 
the orbital plane are the radius vector r and the true anomaly /. In terms of r and / one has according to |8q | (see 
also H, chapter 1) 

r = r (P cos / + Q sin /) , (121) 

where the unit vectors P, Q are defined by 

P = lcosw + msino; , Q = — lsinw + mcosu . (122) 

The coordinate velocity of the pulsar's companion is given by 
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m„ . 
— - r , 

M ' 

' CM\ I//2 
. = I — J [-Psin/ + Q(cos/ + e)] , 



(123) 
(124) 



where M = m p + m c ; p = a^(l — e 2 ) 1 ' 2 is the focal parameter of the elliptic orbit, and m p and m c are the masses 
of the pulsar and its companion. Accounting for relationships 



r = o,r(1 — ecosu) , rcos/ = or(cosu — e) , rsinf = or(1 — e ) ' sinu 



(125) 



where u is the eccentric anomaly relating to the time of emission, to = T, and the moment of the first passage of the 
pulsar through the periastron, Tq, by the Kepler transcendental equation 



u — e sin u = nb(T — Tq) 



we obtain 



k • r = — o,r sini 
r ■ v c — ~a c an nj, e sin u , 
k ■ v c — a c ri\, (1 — e )~ ' sini 



(cosu — e) sinuj + (1 — e 2 ) 1 ' 2 coswsinw, 

(cosu — e) coscj — (1 — e 2 ) 1 / 2 sinwsinw, 



e cos cj 



(126) 

(127) 
(128) 
(129) 



1 — e cos u 

Here a c = a^ m p /M, and n& = {GM/a^) 1 / 2 is the orbital frequency related to the orbital period Pb by the equation 
n b = 2ir/P b . 



Ignoring all constant factors, the set of equations given in this section allows to write down the Shapiro delay (IIS 
in the form 
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As(T) = — ln^ 1 — ecosu — sini 



sino;(cosM — e) + (1 — e 2 ) 1 ' 2 coswsinu 



(130) 



2tt x m 



2tt sin i x m 



sin i Ph m 



e sin u — 



b mc 



(l-e 2 y/ 2 P b 



smw 



(cos-u- e) + (1 - e 2 ) 1/2 coswf 



e cos lu 



(cosu — e) cosuj — (1 — e 2 ) 1 ' 2 sinwsinw 



1 — e cos u 

where in front of the logarithmic function we have omitted the term of order (Gm/c 3 )(v c /c) which is small and hardly 
be detectable in future. The term of order X/r in the argument of the logarithm is also too small and is omitted. The 
magnitude of the velocity-dependent terms in the argument of the logarithm is of order 1CP 3 -r- 10~ 4 . These terms can 
be comparable with the main terms in the argument of the logarithm when the pulsar is near the superior conjunction 
with the companion and the orbit is nearly edge-on. The velocity-dependent terms cause a small surplus distortion in 
the shape of the Shapiro effect which may be measurable in future timing observations when better precision and time 
resolution will be achieved. Unfortunately, existence of the, so-called, bending time delay |p3| may make observation 
of the velocity-dependent terms in the Shapiro time delay a rather hard problem. 
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B. Moving Gravitational Lenses 

The theoretical study of astrophysical phenomena caused by a moving gravitational lens certainly deserves a fixed 
attention. Though effects produced by the motion of the lens are difficult to measure, they can give us an additional 
valuable information on the lens parameters. In particular, a lensing object moving across the line of sight should 
cause a red-shift difference between multiple images of a background object like a quasar lensed by a galaxy, and a 
brightness anisotropy in the microwave background radiation p2j . Moreover, velocity-dependent terms in the equation 
of gravitational lens along with proper motion of the deflector can distort the shape and the amplitude of magnification 
curve observed in a microlensing event. Slowly moving gravitational lenses are 'conventional' astrophysical objects and 
effects caused by their motion are small and hardly detectable. However, a cosmic string, for example, may produce 
a noticeable observable effect if it has sufficient mass per unit length. Gradually increasing precision of spectral and 
photometric astronomical observations will make it possible to measure all these and other effects in a foreseeable 
future. 

1. Gravitational Lens Equation 

In this section we derive the equation of a moving gravitational lens for the case that the velocity v a of the a-th 
light-ray-deflecting mass is constant but without any other restrictions on its magnitude. This assumption simplifies 
calculations of all required integrals allowing to bring them to a manageable form. In what follows it is convenient to 
introduce two vectors <; a = x(s) — x a (s) and £oa = x(2to — so) — x a (so) (see Fig. for more details on the geometry 
of lens). We also shall suppose that the length of vector s a is small compared to any of the distances: R = |x — Xo|, 
r a = |x — x a (s)|, or r Qa = |x — x a (s )|. It is not difficult to prove by straightforward calculations, taking account of 
the light-cone equation, that 

S a = r a - kr a , <; 0a = r 0a + kr 0a , (131) 

where, as in the other parts of the present paper, we have r Q = x — x Q (s) and ro a = xo — x a (so)- From these equalities 
it follows that 

"■'• = "£• k -«"&' (132) 

and 

d 2 
r Qa - k • r 0Q = 2r 0a - tt^ > ( 133 ) 

2r 0a 

where distances d a — \<; a \ and d 0a = |soa| are Euclidean lengths of corresponding vectors. We can see as well that 



making use of the relationships (131) yields 
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r a - v a • r a = r a (l - k • v a ) - s a ■ v a = r a (l - k • v a ) + O (v a d a ) , (134) 

and the residual term can be neglected because of its smallness compared to the first one. 

It is worth noting that the vector <; a is approximately equal to the impact parameter of the light ray trajectory 
with respect to the position of the deflector at the retarded time s. Indeed, let us introduce the vectors £ 8 = P^a? 
and £* = P l ,x J a (s) which are lying in the plane being orthogonal to the unperturbed trajectory of light ray. Then, 



from the definitions (131), (132) one immediately derives the exact relationship 



^-^ a = ^ a +k^, (135) 



from which follows 



s a = € - L - k |s- , (136) 

and the similar relationships may be derived for <r 0a - It is worthwhile to note that 

P]rl = P)<? = C - |« , (137) 

and 

d 2 \£-£ I 2 d 4 
ra-k.r = ^ = ^-i^ + |%. (138) 

2r a 2r a 8r£ 

Let us denote the total angle of light deflection caused by the a-th body as 

<(r) = 4mJ-^J -fcj^ . (139) 



Thus, for the vectors a 1 and ff introduced in (|68|), ( Ba ) and from the formulas (|74|)-(|7q) one obtains |9 

/ Gm a Va \ / Gm a d a \ 

V c 2 r a c / V c2r a r- a J 



*\t)=J2<{t) + o(^ V ^) + o(^ I ±) , (140) 



(141) 



«™)=-i|:^M^, + o(^) + o(^), 



where (by definition) the transverse velocity v aT = -P 1 ,-^ is the projection of the velocity of the a-th body onto the 
plane being orthogonal to the unperturbed light trajectory. 

Let us introduce the new operator of projection onto the plane which is orthogonal to the vector K 
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V i] = 5 ij _ K i K j _ ( 144 ) 

It is worth emphasizing that the operator V^ differs from P 2 - 7 = <5 y — k l k^ by relativistic corrections because of the 
relation (|37]) between the vectors k and K. We define a new impact parameter £ l = V)x? = VI- Xq of the unperturbed 
light trajectory with respect to the direction defined by the vector K. The old impact parameter £ differs from the 
new one £ by relativistic corrections. The direction of the perturbed light trajectory at the point of observation is 
determined by the unit vector s according to equation (164) . We use that definition to draw a straight line originating 
from the point of observation and directed along the vector s up to the point of its intersection with the lens plane 
(see Fig [3|) . The line is parametrized through the parameter A and its equation is given by 

x l (X)^x l (t) + s i (X-t) , (145) 

where A should be understood as the running parameter, t is the value of the parameter A fixed at the moment of 
observation, and x l (t) are the spatial coordinates of the point of observation. On the other hand, the coordinates of 



the point x l (X) at the instant of time A* when the line (145) intersects the lens plane, can be defined also as 



x i (X*)=X i (X*)+r ) i -&, (146) 

where rf = V\ a; 1 (A*) is the perturbed value of the impact parameter £ l caused by the influence of the combined 
gravitational fields of the (micro) lenses m a , X l (X*) = M.^ 1 ^2 a=1 m a x l a (X*) are coordinates of the center of mass of 
the lens at the moment A*. When the line ( |145| ) intersects the lens plane the numerical value of A up to corrections 
of order 0(d/r) is equal to that of the retarded time s defined by equation like (UAh in which r a is replaced by r - 
the distance from observer to the lens. It means that at the lens plane A* — t ~ — r. Accounting for this note, and 
applying the operator of projection Vij to the equation ( |145| ), we obtain 

rf = e - [a\r) + ^(r) - P\r ) + 7 *(r)] r . (147) 



Finally, making use of the relationships (14C)-(142) and expanding distances r a , r$ a around the values r, r$ respectively 
(see Fig. |5| for explanation of meaning of these distances), the equation of gravitational lens in vectorial notations 
reads as follows 

„ = $-^°a(0 + £«(€), (148) 



where 



A ; 



i-k- Va e-e a 



" (e -'£"-7Ptif^ 
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(149) 
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m * V aT , (\£-£a 



fri x/l - vl \ 2r a r 0a J 
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It is not difficult to realize that the third term on the right hand side of the equation (148) is (d a /ro)(v a /c) times 
smaller than the second one. For this reason we are allowed to neglect it and represent the equation of gravitational 
lensing in its conventional form [M , |9qj 

r? = £-^«(0, (151) 



where ct(£) is given by (149). It is worthwhile emphasizing that although the assumption of constant velocities 



of particles v a was made, the equation (151) is actually valid for arbitrary velocities under the condition that the 



accelerations of the bodies are small and can be neglected. 



It is useful to compare the expression for the angle of deflection a 1 given in equation (14E) with that derived 
one in our previous work H. In that paper we have considered different aspects of astrometric and timing effects 
of gravitational waves from localized sources. The gravitational field of the source was described in terms of static 
monopole, spin dipole, and time-dependent quadrupole moments. Time delay and the angle of light deflection a 1 in 
case of gravitational lensing were obtained in the following form |l| 

t-to = |x-x | -4V> + 2Mm(4rro) , 0^ = 4^, (152) 

where the partial ('projective') derivative reads d{ = P\d/d^ , and r and r$ are distances from the lens to observer 
and the source of light respectively. The quantity ip 1S the, so-called, gravitational lens potential JM| , |9Jj] having the 
form 0j 



fl> = 



M + e im k p S q d j + -I p<i {t*)d pq 



In |C|, (153) 



and €j pq is the fully antisymmetric Levi-Civita symbol. The expression (153) includes the explicit dependence on the 
static mass M, spin S l , and time-dependent quadrupole moment T lJ of the deflector taken at the moment t* of the 
closest approach of the light ray to the origin of the coordinate system which was chosen at the center of mass of 
the deflector emitting gravitational waves so that the dipole moment X % of the system equals to zero identically. It 
generalizes the result obtained independently in |ll]] for the case of a stationary gravitational field of the deflector for 
the gravitational lens potential which is a function of time. In case of the isolated astronomical system of N bodies 
the multipolc moments arc defined in the Newtonian approximation as follows 

JV N N N , v 

M = Yl ma ' ^ = Y1 maX * ' S * = X! m °( Xa X Vq )* ' X V = X! m ° [faXi - 3 X a s ij J ) (154) 

a— 1 a— 1 a— 1 a— 1 ^ ^ 

where the symbol 'x' denotes the usual Euclidean cross product and, what is more important, coordinates and 
velocities of all bodies are taken at one and the same instant of time. In the rest of this section we assume that 
velocity of light-ray-deflecting bodies are small and the origin of coordinate frame is chosen at the barycenter of the 
gravitational lens system. It means that 

N N 

T{t)=Y^rn a x i a {t) = Q, and V (t) = ^m a v l a (t) = . (155) 



a=l 
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Now it is worthwhile to note that coordinates of gravitating bodies in ( J149| ) are taken at different instants of retarded 
time defined for each body by the equation ( Jll| ) . In the case of gravitational lensing all these retarded times are close 
to the moment of the closest approach t* and we are allowed to use the Taylor expansion of the quantity 



N N N 



£ m a xi(s) = J2 m a xi(t*) + J2 m a vi(t*)(s -t*) + 0(s - t* f . (156) 

a— 1 a— 1 a— 1 

Remembering that retarded time s is defined by equation (O) and the moment of the closest approach is given by 
the relationship 

t* = t-k-x = f-k-r -k-Xa(a) , (157) 



we obtain, accounting for (^38|) 



2r a 
Finally, we conclude that 



a - t* = k • x a (s) - £- ~ k- x Q (r) + O (£) + O (^ x a ) . (158) 



N N 



J2maK(s) = J2 mavUnik ■ x Q (i*)] + ... , (159) 



where ellipses denote terms of higher order of magnitude, and where the equation (155) has been used. 

Let us assume that the impact parameter £* is always larger than the distance £* . Then making use of the Taylor 
expansion of the right hand side of equation ( |149| ) with respect to £* and v a /c one can prove that the deflection angle 
a 1 is represented in the form 

on - 49,* , (160) 

where the potential ^f is given as follows 

( N N N 

* = < ^2 m a -k-^m Q v a (s) -y^m a xi(s)dj + (161) 

^a— 1 a— 1 a— 1 

N _, N 



1 1 

k-^m a v a (s)4(s) d 3 + -^2m a x p a (s) x q a (s) d pq Un\£\ 

a=l a=l ' 



and ellipses again denote residual terms of higher order of magnitude. Expanding all terms depending on retarded 
time in this formula with respect to the time £*, noting that the second 'projective' derivative d pq is traceless, and 



taking into account the relationship (159), the center-of-mass conditions (155), the definitions of multipole moments 
( |154| ), and the vector equality 

4(k ■ v Q ) - vi(k ■ x Q ) = (k x (x a x v a )y , (162) 

we find out that with necessary accuracy the gravitational lens potential is given by ^ = tp J56[ . Hence, the gravita- 
tional lens formalism elaborated in this paper gives the same result for the angle of deflection of light as it is shown 



in formulas (|152|), ( |153| ). 
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2. Gravitational Shift of Frequency by a Moving Gravitational Lens 



We assume that the velocity v a of each body comprising the lens is almost constant so that we can neglect the 
bodies' acceleration as it was assumed in the previous section. The calculation of the gravitational shift of frequency 
by a moving gravitational lens is performed by making use of a general equation (fq) . As we are primarily interested in 
gravitational lensing, derivatives of proper times of the source of light, 7o, and observer, T, with respect to coordinate 
time, t, can be calculated neglecting contributions from the metric tensor. It yields 



dt 



= Vl-«5> 



(163) 



dt 



Accounting for the identity (p7[), 

1 + K • v - 2^ m a 



dT y/T~vi 
obtain from (|8l|) 



;v 



dtp 

dt 



ds d dt* d dk % d 
lhlh + ~dtdt* + dt dk { 



B a (s.s 0l t* 1 'k) 



N 



1 + K • v + 2^ m a 



ds d ds d dt* d dk l d 
~dT ~d~s + dt ds + ~dT ~dt* + dt d¥ 



-B Q (s,s ,t*,k) 



After taking partial derivatives with the help of relationships (|84|)-(|93J), using the expansions ( |13 
neglecting terms of order d a /r a , m a /r a , m a /raa, and reducing similar terms, one gets 

1 - k • v a (k x v) • (k x r a ) 



l_„2xl/2 



1 + (K + [3 - p ) ■ v + 4V m a -. _ 

~, \/l - v 



\t-£a\ 



N 



l + (K + /3-/3 )-v +4^ m a —== 

a =i V 1 - v - 



1 - k • v (k x v a ) • (k x r a ) 



\i-U 
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(134), (138) 



(166) 



where the relativistic corrections (3 = /9(r, £), /3 = /3(tq,£) are given by means of expressions (|38|), (|39|), ( |70|) - (|7^). 
Making use of relationship (p7|) between the unit vectors K and k, the previous formula can be displayed as follows 

, , , A ST 1 - k ■ v a (k x v) • (k x r a ) 

X1 2x i/ 2 l-k-v + 4^ m a 

1 + z - 



1-v 2 
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k- v 



r-> 1 — K ■ V 

4 <L m ° n : 

r~f \/l-U; 



1 - k ■ v (k x v a ) • (k x r a ) 



(167) 



This formula is gauge-invariant with respect to small coordinate transformations in the first post-Minkowskian ap- 



proximation which leave the coordinates asymptotically Minkowskian. Moreover, the formula (167) is invariant with 
respect to Lorentz transformations and can be applied for arbitrary large velocities of observer, source of light, and 



gravitational lens. In case of slow motion of the source of light, the equation (166) can be further simplifed by expan- 
sion with respect to powers of Vo/c, v/c, and v a /c. Neglecting terms of order w 4 /c 4 , Vq/c 4 , (m a /d a )(v 2 /c 2 ), etc., this 
yields for the frequency shift 
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(168) 



,\r^ 1 -k- v a [kx (v- v a )] • (k x r a ) 
4 2^ m ° n 2 



where 5v = v — vq. The terms on the right hand side of this formula depending only on the velocities of the source 
of light and observer are the part of the special relativistic Doppler shift of frequency caused by the motion of the 



observer and the source of light. The last term on the right hand side of (168) describes the gravitational shift 
of frequency caused by the time-dependent gravitational deflection of light rays due to relative motion of lens with 



respect to observer [100 . It shows that the static gravitational lens being at rest with respect to observer does not 
lead to the gravitational shift of frequency which appears only if there is a relative transverse velocity of the lens 
with respect to the observer which brings in the dependence of the impact parameter for the a-th body, £ — £ a , on 



time 1 101 . By expanding the last term in the expression (168) with respect to powers £/£ a the gravitational shift of 
frequency reads 



5v\ dip ... 

^)_= 4 5F +V ' a(€) 



(169) 



y 



where the deflection angle a is displayed in (152). It is remarkable that the formula (169) is a direct consequence of 



the equation (152) for the time delay in gravitational lensing. 



Indeed, let us assume that the lens is comprised of an ensemble of N point-like bodies each moving with (time- 
dependent) velocity v a with respect to the origin of the coordinate system chosen near the (moving) barycenter of 
the lensing object. The velocity V of the center-of-mass of the gravitational lens is defined as the first time derivative 



of the dipole moment X % of the lens shown in (154), that is 



M 



a=l 



V 
M 



(170) 



Using this definition and assuming that at the initial epoch the barycenter of the lens is at the origin of the coordinate 



system, we find out that the time derivatives of the lens gravitational potential (28C) reads as follows (sec the remarks 

0d> 



in 1 101 1, [102] which clarify calculation of the derivatives) 



chp_ 
dt* 
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dt 



V l) d i^ = t( v 0-C*) 



(171) 



Formulas given in (171) allow to find out the total differential of equation (152) in terms of the increments of time 



dt and dto . While finding the total differential of the gravitational lens equation, it should be also kept in mind that 
asymptotically in the limit tq — > +oo, r = const., the following relationship (see formulas (p7|), (141), ( 142 )) 



K 



(172) 



between vectors K and k holds. Taking differential of equation (152) using results of (171) and (172) one confirms 



the validity of the presentation (16E) for the gravitational shift of frequency in gravitational lensing. Formula (11691) 
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reflects the fact that the gravitational shift of frequency can be induced if and only if the gravitational lens potential 
is a function of time. 

One sees that, in general, not only the translational motion of the lens with respect to observer generates the 
gravitational shift of frequency but also the time-dependent part of the quadrupole moment of the lens. Besides 
this, we emphasize that the motion of observer with respect to the solar system barycenter should produce periodic 
annual changes in the observed spectra of images of background sources in cosmological gravitational lenses. This 



is because of the presence of the solar system time-periodic part of velocity v of observer in equation (169). The 



effect of the frequency shift may reveal the small scale variations of the temperature of the CMB radiation in the sky 
caused by the time-dependent gravitational lens effect on clusters of galaxies having peculiar motion with respect to 
the cosmological expansion. However, it will be technically challenging to observe this effect because of its smallness. 



The simple relationship (169) can be compared to the result of the calculations by Birkinshaw & Gull ( [103|, 
equation 9) . We have checked that the derivation of the corresponding formula for the gravitational shift of frequency 



given by Birkinshaw & Gull [ 103 1 on the ground of a pure phenomenological approach and cited in |92] is consistent, 



at least, in the first order approximation with respect to the velocity of lens [104]. Preliminary numerical simulations 



of the CMB anisotropies by moving gravitational lenses carried out in the paper 106] on the premise of formula (1169) 
under assumption v = 0, confirm a significance of the effect for future space experiments being designed for detection 
of the small scale temperature fluctuations of the CMB. 

However, we would like to make it clear that in practice the gravitational shift of frequency caused by moving 



gravitational lens must be calculated on the basis of a different-from-equation (169) formulation. The matter is that 



the gravitational lens is not at infinity but at a finite distance. For this reason, calculation and subtraction of special 
relativistic Doppler shift of frequency in equation ( |16§| ) should be done using the unit vector K related to k by the 
transformation (p7[). Using the given transformation for replacement of k by K in ( 169 ), remembering (see equations 
( |l4l| ), ( |142| ) ) that /3(t) = —r/R a., and the angle /3(tq) is negligibly small, we obtain for the observable shift of 
frequency in gravitational lensing 



obs 



= 4 ^ + i? (V ' a)+ i? (V0 ' a) ' 



(173) 



where r and rg are distances from observer to the lens and from the lens to the source of light respectively, R 
r + tq. In the limit tq — > +oo, r = const., the equation (U73) goes over to equation (169). 



x-xo 



The last two terms in the right hand side of equation (173) have been derived by Bertotti & Giampieri [107] who 
used a different mathematical technique assuming that the lens is static. Hence, they missed the first term in the right 



hand side of (|173|) discovered by Birkinshaw & Gull |103| who, in their own turn, neglected the contributions due to 



the motion of source of light and observer. It is also useful to note that equation (22) in the paper [107] for Doppler 
shift in gravitational lensing contains a misprint of algebraic sign in front of the term depending on the velocity of 



observer. The error has been corrected in the paper [108] by less et al. (see equation (8) in [ 108 1 ) so that their result 



37 



coincides precisely with the last two terms in the right hand side of our equation (173) 



C. General Relativistic Astrometry in the Solar System 

1. Theoretical Background 

For a long time the basic theoretical principles of general relativistic astrometry in the solar system were based on 
using the post-Newtoninan approximate solution of the Einstein field equations M - M, [p3|, [109]. The metric tensor 



of the post-Newtonian solution is an instantaneous function of coordinate time t. It depends on the field point, x, the 
coordinates, x a (£), and velocities, v (t), of the gravitating bodies and is valid only inside the near zone of the solar 
system because of the expansion of retarded integrals with respect to the small parameter v a /c B. This expansion 
restricts the domain of validity for which the propagation of light rays can be considered from the mathematical point 
of view in a self-consistent manner by the boundary of the near zone. Finding a solution of the equations of light 
propagation (|19| ) in the near zone of the, for instance, solar system can be achieved by means of expanding positions 
and velocities of the solar system bodies in Taylor series around some fixed instant of time, their substitution into 
the equations of motion of photons (|l9|), and their subsequent integration with respect to time. Such an approach 



is theoretically well justified for a proper description of radar HC] and lunar laser |111| ranging experiments, and 



the interpretation of the Doppler tracking of satellites 1 107 , [ 108 1 , [112] - [114|. The only problem which arises in 



the approach under discussion is how to determine that fiducial instant of time to which coordinates and velocities of 
gravitating bodies should be anchored. Actually, the answer on this question is dimmed if one works in the framework 
of the post-Newtonian approximation scheme which disguises the hyperbolic character of the Einstein equations for 
the gravitational field and does not admit us to distinguish between advanced and retarded solutions of the field 
equations pQ]. For this reason, propagation of light rays, which always takes place along the isotropic characteristics 
of a light cone, is different in the post-Newtonian scheme from the gravitational-field propagation because the latter 
propagates in that framework instantcneously and with infinite speed. Thus, the true causal relationship between 
the position of the light particle and location of the light-ray-deflecting bodies in the system is violated which leads 
to a necessity to use some artificial assumptions about the initial values of positions and velocities of the bodies for 
integration of equations of light propagation (see Fig. ^ for more details). One of the reasonable choices is to fix 
coordinates and velocities of the body at the moment of the closest approach of light ray to it. Such an assumption 
was used by Hellings [ jlOj ] and put on more firm ground by Klioner & Kopeikin [H revealing that it minimizes the 
magnitude of residual terms of the post-Newtonian solution of the equations of light propagation. 

As it has been explained above, the post-Newtonian approach has stringent limitations when applied to the in- 
tegration of equations of light propagation in the case when the light-ray-perturbing gravitating system is not in a 
steady state and the points of emission, Xo, and observation, x, of light are separated by the distance which is much 
larger than the characteristic (Keplerian) time of the system. The first limitation comes from the fact that the general 
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post-Newtonian expansion of the metric tensor diverges as the distance r from the system increases (see, for instance, 



[115 1 - [117]). Usually this fact has been ignored by previous researches who used for the integration of the equations 
of light rays the following trancated form of the metric tensor 

1 , 2C/(i,x) 



000 (*,x) 



+ 0(c" 4 ) 



(174) 



0oi(*.x) 



4£/ 4 (t,x) 
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where the instantaneous, Newtonian-like potentials are given by the expressions 
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(176) 



(177) 



(178) 



and all terms describing high-order multipoles have been omitted [118]. From a purely formal point of view the 



expressions ( |174[ )-(176) are not divergent when the distance r approaches infinity but residual terms in the metric 
tensor are. It means that the post-Newtonian metric can not be used for finding solutions of equations of light 
propagation if the distance r is larger some specific value tq. This spatial divergency of the metric tensor relates to 



the fact that the expressions (174)-(|l7q) represent only the first terms in the near-zone expansion of the metric and 



say nothing about the behavior of the metric in the far-zone Iq], [117|, [119]. The distance ro bounding the near-zone 
is about the characteristic wavelength \ gr of the gravitational radiation emitted by the system (A gr ~ (cPb )/(47r), 
where P& is the characteristic Keplerian time of the system). If we assume, for example, that the main bulk of the 
gravitational radiation emitted by the solar system is produced by the orbital motion of Jupiter, the distance ro does 
not exceed 0.3 pc. Almost all extra-solar luminous objects visible in the sky lie far beyond this distance. From this 
point of view, the results of integration of the equations of light propagation from stars of our galaxy and extra-galactic 
objects having been performed previously by different authors on the premise of the implementation of metric tensor 



(174)-( |l76| ) can not be considered as rigorous and conclusive for residual terms of such an integration were never 
discussed. 

The second limitation for the application of the near-zone expansion of the metric tensor relates to the retarded 
character of the propagation of the gravitational interaction. The expressions (|l74[)-(|l76|) are instantaneous functions 



of time and do not show this property of retardation at all. At the same time the post- Newtonian metric (174) 



(176) can be still used for integration of equations of light rays, at least from the formal point of view, because the 
integration will give a convergent result. However, we may expect that the trajectory of light ray obtained by solving 



:-!!) 



the equations of propagation of light using the instantaneous potentials will deviate from that obtained using the 
metric perturbations expressed as the Lienard- Wiechert potentials. Such a deviation can be, in principle, so large 
that the error might be comparable with the main term of relativistic deflection of light and/or time delay. None of 
the methods of integration attempted so far contains error estimates in a precise mathematical sense; at best, errors 
have been roughly estimated using matched asymptotic technique H . None of other previous authors have ever tried 
to develop a self-consistent approach for calculation of the errors. 

One more problem relates to the method of performing time integration of the instantaneous potentials along the 
unperturbed trajectory of the light ray. This is because coordinates and velocities of bodies are functions of time. 
Even in the case of circular orbits we have a problem of solving integrals of the type 

ft N ft fa 

\ U(t,x)dt=ym a i 7—, s 7—, -, s -, 7TT > (179) 

J t f~[ J t l x o + k (t - t ) - A a |ei am(u a t + ip a ) + e 2 cos{uj a t + ip a )\ \ 

where A a , uj a , and ip a are the radius, the angular frequency, and the initial phase of the orbit of the a-th body 

respectively, and ei, e 2 are the unit orthogonal vectors lying in the orbital plane. The given integral can not be 

performed analytically and requires the application of numerical methods. In case of elliptical motion calculations 



will be even more complicated. Implicitly, it was usually assumed that the main contribution to the integral (17£) 
comes from that part of the trajectory of light ray which passes nearby the body deflecting light rays so that one is 
allowed to fix the position of the body at some instant of time which is close to the moment of the closest approach 
of the light ray to the body. However, errors of such an approximation usually were never disclosed except for the 
attempt made by Klioner & Kopeikin Q|. Nevertheless, it is not obvious so far that the error analysis fullfiled in Q] 
is complete and that the use of the Taylor expansion of coordinates and velocities of the solar system bodies with 
respect to time, made in the neighborhood of the instant of the closest approach of photon to the light-deflecting body 



in order to perform the integration in (179), minimizes errors of calculations and allows to solve light ray equations 



with better precision. Moreover, such an expansion is allowed only if photon moves near or inside of the gravitating 



system. Far outside the system the other method of solving the integral (179) is required [12C j. 

Regarding this difficulty Klioner & Kopeikin H have used a matched asymptotic technique for finding the perturbed 
trajectory of the light ray going to the solar system from a very remote source of light like a pulsar or a quasar. The 
whole space-time was separated in two domains - the near and far zones lying correspondingly inside and outside 
of the distance r$ being approximately equal to the characteristic length of gravitational waves emitted by the solar 
system. The internal solution of the equations of light rays within the near zone have been obtained by expanding 
coordinates and velocities of the bodies in the Taylor time series and then integrating the equations. The external 
solution of the equations has been found by decomposing the metric tensor in gravitational multipoles and accounting 
only for the first monopole term which corresponds mainly to the static, spherically symmetric field of the Sun. A 
global solution was obtained by matching of the internal and external solutions at the buffer region in order to reach 
the required astrometric accuracy of 1 /iarcsec. The approach we have used sounds reasonable and may be used in 
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theoretical calculations. However, it does not help very much to give a final answer to the question at which moment 
of time one has to fix positions and velocities of bodies when integrating equations of light propagation inside the 
near zone. In addition, the approach under consideration does not give any recipe how to integrate equations of light 
propagation in the external domain of space (beyond ro) if the higher, time-dependent gravitational multipolcs should 
be taken into account and what magnitude of the perturbative effects one might expect. In any case, the global 
solution obtained by the matched asymptotic technique consists of two pieces making the visualization of the light 
ray trajectory obscure and the astrometric implementation of the method impractical. 

For these reasons we do not rely in this section upon the technique developed in j| but resort to the method 
of integration of light ray equations based on the usage of the Lienard-Wiechert potentials. This method allows to 
construct a smooth and unique global solution of the light propagation equations from arbitrary distant source of 
light to observer located in the solar system. We are able to handle the integration of the equations more easily and 
can easily estimate the magnitude of all residual terms. Proceeding in this way we also get a unique prediction for 
that moment of time at which coordinates and positions of gravitating bodies should be fixed just as the location of 
observer is known. We shall consider three kinds of observations - pulsar timing, very long baseline intcrferometry 
(VLBI) of quasars, and optical astrometric observations of stars. 

2. Pulsar Timing 

The description of the timing formula is based on the usage of equations ( fi"3| ) and (E3). Taking in the equation 
(p4|), which should be compared with its post-Newtonian analogue ( |179| ), terms up to the order v a /c inclusively, we 
obtain 



B a (s,s Q ) = -In 



r a (s) -k-r a (s) 



k-v Q (CWC 



*t. +k .,.(o-c'- W' (180) 



_r a (so) ~ k ■ r Q (s ) 

where the retarded times s and sq should be calculated from equations flll| ) and ( |46D respectively, r a (s) = x — x a (s), 
r a(so) = xo — x a (so), and we assume that the observation is made at the point with spatial barycentric coordinates x 
at the instant of time t, and the pulsar's pulse is emitted at the moment to from the point xo which is at the distance 
of the pulsar from the solar system, typically more than 100 pc. 

In principle, the first term in this formula is enough to treat the timing data for any pulsar with accuracy required 
for practical purposes. The denominator in the argument of the logarithmic function is r Q (so) — k- r a (so) — 2i?, where 
R is the distance between the barycenter of the solar system and the pulsar. The logarithm of 2R is a function which 
is nearly constant but may have a secular change because of the slow relative motion of the pulsar with respect to the 
solar system. All such terms are absorbed in the pulsar's rotational phase and can not be observed directly. For this 



reason, in what follows, we shall omit the denominator in the logarithmic term of equation (180). Wc emphasize that 



positions of the solar system bodies in the numerator of the logarithmic term are taken at the moment of retarded 
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time which is found by iterations of the equation s = t — |x — x a (s)|. It makes calculation of the Shapiro time delay 
in the solar system theoretically consistent and practically more precise. 

There is a difference between the logarithmic term in (|180|) and the corresponding logarithmic terms in timing 



formulas suggested by Hellings [[L0| and Doroshenko & Kopeikin [ 121 1 where the position of the a-th body is fixed 
at the moment of the closest approach of the pulse to the body. It is, however, be not so important in practice as 
timing observations are not yet precise enough to distinguish the Shapiro delay when positions of bodies are taken at 
the retarded time or at any other one, being close to it. Indeed, the maximal difference is expected to be of order of 
(4GMq/c 3 )(vq/c)(1 + cosfl) -1 , where M Q and v© are mass of the Sun and its barycentric velocity respectively, and 
6 is the angle between directions towards the Sun and the pulsar. For vq is less than 20 m/s and 9 can not exceed 
0.25 degrees the error in timing formula relating to various definitions of the Sun's barycentric coordinates in the 
expression for the Shapiro time delay is less than 200 nanoseconds. Although this value is yet beoynd of observational 



limit it would be desirable to update existing timing data processing programs like TEMPO 122] and TIMAPR [|33| 



[121] to make their functional structure be in agreement with the latest theoretical developments. 



The integral in equation (18C) can not be calculated analytically if the trajectory of motion of the bodies is not 
simple. The matter is that in the case when light propagates from the remote source to the gravitating system the 
time interval between s and so is not small as it was in the case of the derivation of timing formula for binary pulsars 
in section 7A (see equation ( p6| ) and related discussion as well as caption to Fig. 0) . This was because light propagates 
from the binary system in the same direction as gravitational waves emitted by it, so that the gravitational field of the 
binary system is almost "frozen" as seen by the outgoing photon. When we consider propagation of light towards the 
solar system the infalling photon moves in the direction being opposite to that of propagation of gravitational field 
generated by the moving solar system bodies. For this reason, the difference s — sq ~ 2R, and is very large. Thus, 
we are not allowed, as it was in the case of derivation of timing formula for binary pulsars, to use the expansion of 
coordinates and velocities of the solar system bodies in Taylor series with respect to time. Moreover, integrals, like 



that in (18C), should be also calculated without any expansion using the known law of motion of gravitating bodies, 
that is the solar system ephemerides like DE200, DE245, or an equivalent one. Let us give an idea what kind of result 
we can get proceeding in this way. 

First of all, we note that the orbital plane of any of the solar system bodies lies very close to the ecliptic and can 
be approximated fairly good by circular motion up to the first order correction with respect to the orbital eccentricity 
which is usually small. The motion of the Sun with respect to the barycenter of the solar system may be described 
as a sum of harmonics corresponding to gravitational perturbations from Jupiter, Saturn, and other smaller bodies. 
Thus, we assume that x Q is given in the ecliptic plane as follows 

x a (i) = A [cos(ni) ei + sin(nt) e 2 ] , (181) 

where A and n are the amplitude and frequency of the corresponding harmonic in the Fourier decomposition of the 
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orbital motion of the a-th body, ei is directed to the point of the vernal equinox, e 2 is orthogonal to ei and lies in 
the ecliptic plane. The vector k is defined in ecliptic coordinates as 

k = — cos b cos / ei — cos b sin I e 2 — sin b e 3 , (182) 

where b and I are ecliptic spherical coordinates of the pulsar. Substituting these definitions into the integral of equation 
( |180| ) and performing calculations using approximate relationship £ = t* — y + k • x. a (y), where y is the new variable 
defined in (f48|), we get 

(' +*l'J ai0 (^ r = k ■ v «(**)( Ci [« (^ - k • r B )] - Ci [n (r 0a - k • r 0a )] 1 (183) 

-k-x a (i*K Si [n(r a - k • r )] - Si [n (r 0a -k-r 0o )j 

where t* = t — k • x is the time of the closest approach of the light ray to the barycenter of the solar system. Taking 

into account the asymptotic behaviour of sine and cosine integrals for large and small values of their arguments in 

relationship (Q83) we arrive at the approximate formula for the Shapiro delay 

N (C \ 

A(«) = -53m„[l-k.v a (O]ln[ra(*)-k.r a (s)] + 0f^^J , (184) 

a— 1 ^ ' 

where the residual term denotes all contributions which are simple products of the gravitational radius Gm a /c 3 of 
the a-th body, expressed in time units, by the ratio v a /c up to a constant factor. If one takes numerical values of 
masses and velocities of the solar system bodies one finds that such residual terms are extremely much smaller than 
the level of errors in timing measurements. We conclude that these residual terms can not be detected by the present 
day pulsar timing techniques. 

3. Very Long Baseline Interferometry 

VLBI measures the time differences in the arrival of microwave signals from extragalactic radio sources received at 
two or more radio observatories ]123[ |. Generally, geodetic observing sessions run for 24 hours and observe a number 
of different radio sources distributed across the sky. The observatories can be widely separated; the sensitivity of 
the observations to variations in the orientation of the Earth increases with the size of the VLBI network. VLBI 
is the only technique capable of measuring all components of the Earth's orientation accurately and simultaneously. 
Currently, VLBI determinations of Earth-rotation variations, and of the coordinates of terrestrial sites and celestial 



objects are made routinely and regularly with estimated accuracies of about +/-0.2 milliarcsecond or better [ 123 



[124|. Such a high precision of observations requires an extremely accurate accounting for different physical effects in 
propagation of light from radio sources to observer including relativistic gravitational time delay. 

There have been many papers dealing with relativistic effects which must be accounted for in VLBI data processing 



(see, e.g., [109|, [ 125 1 , [126], and references therein). The common efforts of many researches in this area have resulted 
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in the creation of what is commonly believed now to be as a 'standard' model of VLBI data processing which is called 



a consensus model [126 emerged from a workshop held in 1990 |12q [. The accuracy limit chosen for the consensus 
VLBI relativistic time delay model is 10 -12 seconds (one picosecond) of differential VLBI delay for baselines less than 
two Earth radii in length. As it was stated, in the model all terms of order 10 -13 seconds or larger were included 
to ensure that the final result was accurate at the picosecond level. By definition, extragalactic source coordinates 
derived from the consensus model should have no apparent motions due to solar system relativistic effects at the 
picosecond level. Our purpose in this section is to analyze critically this statement and to show that the consensus 
model is not enough elaborated, at least theoretically, in accounting for relativistic effects in propagation of light at 
the picosecond level. For this reason, we propose necessary modification of the consensus model to make it applicable 
at the level of accuracy approaching to 10~ 13 seconds without any restrictions. 

In what follows we work for simplicity with the barycentric coordinate time of the solar system only. Precise 
definition of the measuring procedure applied in VLBI requires, however, derivation of relativistic relationship between 



the proper time of observer and the barycentric coordinate time. It is given, e.g. in [127], and can be added to the 
formalism of the present section for adapting it to practical applications. A complete description of such an extended 
formalism will be given elsewhere. 

The VLBI time delay (see Fig. ||) to be calculated is the time of arrival of electromagnetic signal, t 2 , at station 
2 minus the time of arrival of the same signal, £1, at station 1. The time of arrival at station 1 serves as the time 
reference for the measurement. In what follows, unless explicitly stated otherwise, all vectors and scalar quantities are 
assumed to be calculated at t\ except for position of the source of light xo which is always calculated at the of time 
of light emission to- We use for calculation of the VLBI time delay equations ([43]), (J44|) referred to the barycentric 
coordinate frame of the solar system. The equations give us 

t 2 -ti = |x 2 (t 2 ) - x | - |xi - xo| + A(t 2) *o) - A(ti, t ) , (185) 

where Xo are coordinates of the source of light, x 2 (£ 2 ) are coordinates of the station 2 at the moment t 2 , Xi are 
coordinates of the station 1 at the moment t\ . The differential relativistic time delay is given in the form 

N 

A(t 2 ,t ) - Ait^to) = 2j2[B a (s2,s ) - B a { Sl ,s )} , (186) 

0=1 

where the difference of the B a 's up to the linear with respect to velocities of the solar system bodies reads (see equation 

B a (s 2 , s ) - B a ( Sl , so) = In ria ~ ^ " ria - In r ° a " ^ ' T ° a (187) 

r 2a - k 2 • v 2a r 0a - k 2 • r 0a 

+k 2 • v a (s 2 ) In (r 2a - k 2 • r 2a ) - k : • v a (si) In (r la - k x ■ r lo ) 
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+ r In [t$ + k! • x Q (C) - C] ki • v o (0 d( ~ r In [i* + k 2 • x a (C) - C] k 2 • v„(C) d( • 

Js Js 

Herein s± and s 2 are retarded times determined iteratively from the equations 

si=ti-|x 1 -x (si)|, (188) 

s 2 = *2-|x 2 (i 2 )-x (s 2 )|, (189) 

the quantity ri a = Xi — x Q (si) is the vector from the a-th body to the station 1, r 2a = x 2 (£ 2 ) — x a (s 2 ) is the vector 
from the a-th body to the station 2, r la = |r la |, r 2a = |r 2a | and 

t{ = h - k x • xi , (190) 

^=i 2 -k 2 -x 2 (i 2 ), (191) 

are the moments of the closest approach of the light rays 1 and 2 to the barycenter of the solar system. It will be also 
helpful in comparing our approach with the consensus model to use the moments of the closest approach of the light 
rays 1 and 2 to the a-th body which we will define according to the rule 

*I = ti-k 1 T l0> (192) 

t* 2a = t 2 - k 2 • r 2a . (193) 

It is worth emphasizing that our definitions of times £ la and £ 2a are slightly different from the definitions of similar 
quantities given in the 'standard' consensus model. It relates to the definition of positions of bodies in the vectors 
i"i Q and r 2a . In our case we refer the coordinates of the bodies to the retarded times s\ and s 2 respectively while 
in the consensus model they are taken at the times t\ and £ 2 . It introduces some uncertainty into the notion of the 
instant of the closest approach of light ray to to the body or barycenter of the solar system which appears due to not 
covariant formulation of the relativistic time delay in the consensus model. There is no such an uncertainty in our 
approach which is fully covariant in the first post-Minkowskian approximation. 
The unit vectors ki and k 2 are defined as 

kl = ^V k 2= | X fj- X ° (194) 

|xi-x | |x 2 (i 2 )-x | 

which shows that they have slightly different orientations in space. Let us introduce the barycentric baseline vector 
at the time of arrival t\ through the definition B = x 2 (^) — Xi(ti). Let us stress that the baseline vector lies on the 
hypersurface of constant time £1. The original version of the relativistic relationship of the barycentric baseline vector 



to the geocentric one, b, can be found in [127] or later publications [109|, [126]. We shall neglect this relativistic 
difference in the expression for the Shapiro time delay because it is inessential in our present discussion. Thus, we 
assume B = b. The difference between the vectors ki and k 2 may be found using the expansion with respect to 
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powers of the small parameter b/R where R is the distance between the barycenter of the solar system and source of 
light. We have 

(195) 



x 2 - x = xi - x + b + v 2 (*2 - *i) + O I -s b 



|x 2 - x | = |xi - x | + b ■ ki + v 2 ■ ki (t 2 -h)+0 



O 



R 



where v 2 is the velocity of station 2 with respect to the barycenter of the solar system. These expansions yield 



k 2 = ki 
and for the time delay ( |l85| ) 

t 2 - h = k x • b 



ki x (b x ki) 
R 



O 



1 + v 2 ■ ki + O 



O 



cR 



R 



O 



R 2 



•A(t a ,«o)-A(ti,to) 



As a consequence of the previous expansions we also have the following equalities 

(b x ki)(ki x x 2 ) 



i 2 — i l — 



R 



O 



°'« 



t* -t* 

la la 



(ki • b)(ki ■ v ) - {r 2a - ri a )(ki • v a ) 



(b x ki)(ki x r 2a ) 
R 



O 



O 



R 



(196) 



(197) 



(198) 



(199) 



which evidently shows that, e.g. for the Jupiter and for the source of light at infinity, the time difference i 2Q — i* Q 
is of the order (Rq/c)(vj/c) ~ 75 nanoseconds, that is, rather small but still may be important in the analysis of 
observational errors. For VLBI observations of the solar system objects the time difference t 2Q — i* Q can approach 
the value R®/c ~ 30 ms which can not be ignored at all. The time difference t 2 — t* can be considered for extra-solar 
objects as negligibly small since it is of the order (R^/c) by the annual parallax of the source of light which makes it 
much less than 1 picosecond. In case of VLBI observations of the solar system objects the time difference t 2 — t\ can 



not be ignored anymore but we do not elaborate it here. Now we can simplify formula ( 187 ) 
First of all, taking into account the relationship ( |l96| ), we obtain 

r Qa - ki • r 0a 



In- 



(b x ki)(ki x r 0a ) ^ fb_ 
R(r 0a - ki • r 0a ) ~ \R 



(200) 



r 0a - k 2 • r 0a 
which is of the order of the annual parallax of the source of light. This term can be neglected in the delay formula 



(187) since it gives a contribution to the delay for extra-solar system objects much less than 1 picosecond [128|. After 



noting that in the expression for the difference of the two integrals in (187) one can equate k 2 = ki and t 2 = t*, we 
state that the difference reads 

J 1 In [t\ + kx • x a (C) - C] ki • v a (CR = kx • v Q ( Sl )|(r 2a - k x • r 2a ) ln(r 2a - kj • r 2a ) (201) 

- (r la - ki ■ ri ) ln(n a - ki ■ r la ) + r la - r 2a - ki • r la + ki • r 2a > , 
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and after multiplication by the factor 2Gm a /c 3 is much less than 1 picoarcsecond. Hence, we neglect those two 
integrals from the expression for the VLBI delay A(ti,t2) = A(t2,to) — A(ti,to)- 

Finally, taking into account that ki = — K, up to the corrections of order of the annual parallax, we get for the 
time delay 

N 

m (~\ -l_ XC . IT \ In 

r 2a + K • r 2a 



A(t 1 ,f a ) = 2y)m g (l + K-v a ) ln r ^ + K-fia (2Q2) 



a=l 

where v a = v a («i), r\ a = |ria|, 7" 2q = |r 2a |, and 



r ia =xi(ii) -Xo(si) , r 2a =x 2 (£ 2 ) -x Q (s 2 ) . (203) 



We emphasize that our formula (202) includes the first correction for the velocity of the bodies deflecting light rays. 



Moreover, there is a difference between the definitions of the vectors Y\ a , r 2a in our model (202) and the consensus 



model (see [126], chapter 12, formula (1)). In our case the coordinates of stations Xi, x 2 are taken at the instants ti, 



t 2 respectively, and the coordinates of the light-deflecting bodies are calculated at the retarded times s±, s 2 defined in 



( 188 ) , ( 189 ) which is a direct consequence of our rigorous approach of the integration of equations of light propagation. 



On the other hand, in the consensus model coordinates of stations are taken also at the instants t\, t 2 but coordinates 



of the a-th body are calculated only at the time t\ a defined in ( 192 ) . Strictly speaking, this prescription can be justified 
only for VLBI observations of the distant, extra-solar objects and is only marginally correct for VLBI observations of 
the solar sytem objects. The prescription to obtain the position of the gravitating body at the time of closest approach 
of the ray path to the body was based on an intuitive guess (see, for instance, JlQ] ). Such a guess gives a rather good 
approximation but can not be adopted as a self-consistent theoretical recommendation in doing practically important 
numerical processing of VLBI observations and, especially, in the dedicated experiments specifically designed to test 



gravitational deflection of light in the solar system [12£] 



If we denote by At grav the VLBI delay in the consensus model, as it is described in the IERS Conventions (see 



[126|, formulas (1), (2) of chapter 12, or formula (5) from [109]), and put the PPN parameter 7 = 1, we get in the 
framework of General Relativity the following relationship between the Lorentz-covariant expression for the time delay 
in our model and At gra v- 

\n *\ a + 1 nV Gm « (v m ria ± K ' ria o V^ Gm a (v„ x r ia )(b x r ia ) 

A(t X , t 2 ) = Mgrav + 2 ^ — — (K • v a ) In - 2 ^ —4 -3 + ■ ■■ , (204) 

a =l c r 2a +i^ r 2a a=i c r la 

where ellipses denote residual terms, and we have restored the universal gravitational constant G and speed of light c 

for convenience. One can see that the Shapiro time delay in the consensus model was not properly defined although 

it had no consequences for practical observations in the recent past. Indeed, the third term in the right hand side of 



( 204 ) is so small that can be neglected for any observational configuration of the source of light and the deflecting 



body including the Earth. Expansion of the second term in the right hand side of equation (204) with respect to 



powers b/d a , where d a is the impact parameter of the light ray with respect to the a-th light-deflecting body, gives 
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2 ^(K.v Q )ln ^ + g- ri " =-2^(K.v a ) b - (ni " + K) =-4^(K.v a ) b - (n " + K) ^, (205) 
c 4 r 2a + K • r 2a c 4 r la + K • r la c 4 d a d a 

where the unit vector ni a = ri /n . For the light ray grazing, for example, the limb of the Sun the term under 
consideration can reach a few picoseconds. The effect amounts 1 picosecond for radio source being at the angular 
distance 10 arcminutes from the Jupiter if Jupiter is at the distance 5 AU from the Earth. This may will have real 
impact in near future on the treatment of the gravitational deflection of light by massive solar system plantes in the 
specialized high-precision VLBI experiments. 

We would like to note that relativistic perturbations of light propagation caused by velocities of moving gravitating 



bodies were considered by Klioner [130] in order to find corresponding corrections to the consensus model of VLBI 
data processing. That author approached the problem doing calculations on the base of the post- Newtonian metric 
tensor. As we have shown in the present paper, such an approximation is not exact enough to take properly into 
account all effects of retardation in the metric which contribute to velocity dependent terms in the propagation of 



light. Nevertheless, at least formally, the result published in |130| , equation (4.9), coincides with our equation (202) 
but the coordinates and velocities of the light-ray-deflecting bodies are taken at the time of the closest approach of 
photon to the a-th body. This produces errors of the order comparable with the last term shown in the right hand 



side of (204) which are negligibly small. Hence, we conclude that the relativistic model of VLBI data processing 



proposed in [ 130 1 is although theoretically incomplete but practically good enough for applications at the level of 



accuracy about one picosecond for astronomical objects with negligibly small parallaxes. 

4- Relativistic Space Astrometry 

Space astrometry is a new branch of fundamental astrometry. Ground-based telescopes may reach the angular 
resolution not better than 0.01 arcseconds. This limits our ability to create a fundamental inertial system on the sky 



[131 1 with the accuracy required for a much better understanding the laws of translational and rotational motions 
of celestial bodies both inside and outside of the solar system. The epoch of the space astrometry began in 1989 
when the HIPPARCOS satellite was successfully launched by Ariane 4 of the European Space Agency on 8 August 
1989. Despite the failure to put the satellite on the intended geostationary orbit at 36,000 Km from Earth the 



astrometric program has been completely fulfilled [132]. As a result the new astrometric catalogue of all stars up to 
13-th stellar magnitude was obtained. It includes about 120,000 stars and has a precision of around 0.002 arcseconds. 
Unfortunately, such high precision can not be retained longer than 10 years because of errors in determination of 
proper motions of stars. For this reason the second analogous mission having the same or better astrometric accuracy 
should be launched in near future. 

Rapid industrial development of space technologies allows us to hope that in the next several years the precision 
of astrometric satellites will reach a few microarcseconds or even better in the determination of positions, proper 
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motions, and parallaxes of celestial objects. All together, the photometric sensitivity of measuring devices will be 
substantially improved. As an example, we refer to a new space project of the European Space Agency named GAIA 
(Galactic Astrometric Interferometer for Astrophysics). In the framework of this project H positions, proper motion, 
and parallaxes of about 1000 million stars up to 20 stellar magnitude arc to be measured with accuracy better than 
10 microarcsecond. It means that practically almost all stars in our Galaxy will be observed and registered. 

Such extremely difficult observations can not be processed adequately if numerous relativistic corrections are not 



taken into account in a proper way [109|. Indeed, the relativistic deflection of light caused by the Sun is not less 



than 1 milliarcsecond throughout all of the sky. Major planets produce a relativistic deflection of light about 1 



microarcsecond at the angular distances from 1 to 90 degrees outside the planet [ 133 1 . It is worth emphasizing that 
the relativistic deflection of light produced by the Earth reaches a maximal value of about 550 microarcseconds and 
should be accounted for any position of a star with respect to the Earth. In addition, the reduction of astrometric 
observations made on the moving platform will require an extremely careful consideration of relativistic aberration 



[134 1 and classic parallax terms in order to reduce the measurements to the solar system barycenter - the point to 
which the origin of the fundamental incrtial system is attached. Perhaps, it would be more properly to say that 
data processing of observations from modern space astrometric satellites should be fully based on general relativistic 
conceptions rather than on a classical approach in which the relativistic corrections are considered as additive and 
are taken into account at the very last stage of the reduction of observations. 

As far as we know, the first attempt to construct such a self-consistent theory of astrometric observations was 



proposed by Brumberg & Kopeikin [135] and further explored in [127] and H. The main idea of the formalism is to 



exploit to a full extent a relativistic theory of reference frames in the solar system developed in papers ( |136| - [139 



rll, and references therein). An independent, but similar approach with more emphasis on mathematical details was 



presented in papers |80j, [140|, [141]. One global and several local reference frames have been constructed by solving 
in a specific way the Einstein equations for gravitational field. The global frame is the barycentric reference frame of 
the solar system with origin at the barycenter. Among the local frames the most important for us is the geocentric 
frame with origin at the geocenter and the proper reference frame of an observer (or the satellite in case of a space 



mission like HIPPARCOS or GAIA). All reference frames are harmonic [142| and were constructed in such a way to 
reduce to minimum all fictitious coordinate perturbations which may be caused by unsophisticated technique in using 
coordinate transformations from one frame to another. We have discovered and outlined corresponding relativistic 
transformations between the frames which generalize the well-known Lorentz transformation in the special theory of 
relativity and minimize the magnitude of unphysical coordinate dependent terms. Proceeding in this way we have 
achieved a significant progress in describing relativistic aberration, classic parallax, and proper motion corrections [H. 
However, the problem of propagation of light rays from distant sources of light to an observer in the non-stationary 
gravitational field of the solar system was not treated thoroughly enough. This section refines the problem and gives 
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its final solution. 

The quantity which we are specifically interested in is the direction towards the source of light (star, quasar) 
measured by a fictitious observer being at rest at the point with the solar system barycentric coordinates (i,x). This 
direction is given, actually, by the equation (f34|) and can explicitly be written as follows 

S -(r,g) = ^ + 2 y^^ (1 'V Va)2 ^ Ti , 2 y^^ 2-k-v a 






j -a 



^VT^2r -k-r a J a R^JTZT^ \ 2R 

a— 1 V a a— 1 V a ^ 

where positions and velocities of the solar system light-deflecting bodies are calculated at the retarded time s = 
t — |x — x |, i? = |x — x a | is the distance from the source of light to observer, and ellipses denote residual terms 
depending on accelerations of the bodies given by the retarded integrals (|74|)-(|76|). We have neglected all terms 
depending on accelerations of the bodies because of their insignificant numerical value. Further simplification of 
equation ( |206|) is possible if we remember that the velocities of bodies, v a , comprising the solar system are small in 
comparison with the speed of light, and distances, R, to stars are very large compared to the size of the solar system. 
This makes it possible to omit all terms being quadratic with respect to v a as well as fifth term in the right hand side 



of (20E) being inversely proportional to R. It yields 

H *\ tsi , ~sr^Gm a ( 2 1 r a \kx(n a xk) r Gm kx(v xk) 
^M^+^-l--k.v a + -v a .n g -- -4^- , (207) 

a— 1 v ' a— 1 

where n a = r a /r a , the sign "x" denotes the usual Euclidean vector product, and we restored the fundamental 



constants G and c for convenience. The equation (207) eliminates incompleteness in the derivation of the similar 
formula given by Klioner p5[ which was obtained using the post-Newtonian expression for the metric tensor and 
under the assumption of rectilinear and uniform motion of the light-deflecting bodies. As we have already noted 
many times in the present paper, the post-Newtonian approximation for the metric tensor does not take into account 



all necessary effects of retardation [143| which are essential in the derivation of the equation (207). Klioner & Kopeikin 
M have simply copied the result of J55| due to the absence at that time of a better theoretical treatment of influence 
of body's velocities on the propagation of light rays. With the mathematical technique invented in the present paper 



the equation (207) gives the correct answer to this question and closes the problem. 



The leading order term in ( 207 ) gives the well-known expression for the angle of deflection of light rays in the 



gravitational field of a static, spherically symmetric body. The velocity dependent terms in ( 207 ) describe small 
corrections which may be important in data analysis of future space missions. The very last term in the large round 
brackets in ( [207| ) may slightly change magnitude of the angle of gravitational deflection for some nearby stars or 
objects within the solar system if the impact parameter of the light ray is small and the deflection angle is expected 
to be rather large. Parallactic corrections to the direction s l are extracted from the unit vector K by its expansion in 
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powers of the ratio (the barycentric distance to observer) /(the barycentric distance to a star). Account for aberrational 
corrections is made by means of relating the direction to the star, s 1 , observed by a fixed fictitious observer, to the 
direction observed by a moving real observer, with the help of the matrix of relativistic transformation displayed in 
section VII of the paper [4]. It is worth emphasizing that the correction for aberration must be done first before 
account for parallax. Complete analysis of the relativistic algorithm of processing observations of celestial objects 
made from a board of a space observatory will be given elsewhere. 

D. Doppler Tracking of Interplanetary Spacecrafts 

1. Approximation Scheme for Calculation of the Doppler shift 



The Doppler tracking of interplanetary spacecrafts [ 144 1 , JL45[, |112|] is the only method presently available to 



search for gravitational waves in the low frequency regime (10 - 1 Hz). Several experiments have been carried 
out so far, for instance, VOYAGER, PIONEER, ULYSSES, GALILEO and MARS-OBSERVER. The space-probe 



CASSINI represents the next step in such gravitational wave Doppler experiments [113|. Its primary target is to 
study the Saturn system. However, the spacecraft carries on board much improved instrumentation and will perform 
three long (40 days each) dedicated data acquisition runs in 2002, 2003 and 2004 to search for gravitational waves 
with expected sensitivity about twenty times better than that achieved so far. The detection of gravitational waves 
requires the precise knowledge of the Doppler frequency shift caused by the solar system's bodies lying near the line 
of sight of observer to spacecraft (see Fig. m. 

Another important implementation of the Doppler tracking is the Global Positioning System (GPS) which uses 
accurate, stable atomic clocks in satellites and on Earth to provide world-wide position and time determination. 
These clocks have relativistic frequency shifts which are so large that, without accounting for numerous relativistic 



effects, the system would not function ( [146], and references therein). Quite recently, the Europeian Space Agency 
(ESA) has adopted a new program aimed at achieving an even better precision in measuring time and frequency in 
space-time observations. The program is called the Atomic Clock Ensemble in Space (ACES) and will be carried out 
on board of the International Space Station (ISS). The principal idea is to use a cold atom clock in absence of gravity 
which will outperform the fountains clock on the ground with the potential accuracy of 5 x 10~ 17 [147]. 



An adequate treatment of such gravitational wave and time-metrology high-precision experiments require advanced 
theoretical development of the corresponding analytic algorithm which properly accounts for all terms of order 10 -16 
and higher in the classic Doppler and gravitational shifts between transmitted and received electromagentic frequencies 
caused by the relative motion of the spacecraft with respect to observer and time-dependent gravitational field of the 
solar system bodies. In this paragraph we discuss basic principles of the Doppler tracking observations and give the 
most important relationships for calculation of the relevant effects. However, the complete theory involves so many 
specific details that it would be unreasonable to give all of them in the present paper. Therefore, only basic elements 
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of the Doppler tracking theory are given here and particular details will be published somewhere else. 

Let us assume (see Fig. |7j) that an electromagnetic signal is being transmitted from the point with barycentric 
coordinates Xo located on the Earth with frequency uq at the barycentric time tp,. It travels to the interplanetary 
spacecraft, is received on its board at the point with barycentric coordinates Xi with frequency V\ at the barycentric 
time t\, and is transponded back to the Earth (on exactly the same frequency v\) where one observes this signal at 
the point with barycentric coordinates X2 with frequency v^ at the barycentric time ti . It is worthwhile to emphasize 
that because of the motion of the receiver with respect to the transmitter during the light travel time of the signal the 
observed frequency V2 is different from the emitted frequency vq even if the signal is transponded from the spacecraft 
being momentarily at rest with respect to the barycentric coordinates of the solar system. 

Proper time of the transmitter at the instant of signal's emission is denoted by 7q and at the instant of the signal's 
reception by 7^. Proper time of the spacecraft's transponder is denoted as T\. Barycentric time at the emission point 
is ioj at the point of reception, £2, and at the specacraft's position, t\. We follow arguments similar to those used in 
section VI. C. The spectral shift of electromagnetic frequency v§ with respect to v\ is given by the equation 

Vn dTi dii dtp, 
v\ dti dtp, dip, 

and the shift of the frequency v\ with respect to V2 is described by the similar relationship 

v-\ dti dto dTo 

V2 dl\ dt\ «E2 

Here the time derivatives dTi/dti and dti/dTi are calculated at the spacecraft's position, dt /dT is calculated at 
the point of emission, and d^jtid at the point of reception. Time derivatives dti/dto and dtijdt\ are obtained from 
the solution of equation of propagation of electromagnetic signal in time-dependent gravitational field of the solar 
system (|43|) establishing theoretical description of the transmitter-spacecraft (up-) and spacecraft-receiver (down-) 
radio links. 

In practice, when Doppler tracking observations are made, the frequency V2 of the receiver is kept fixed. It relates 
to the fact that the frequency band of the receiver must be rather narrow to decrease the level of stochastic noise 
fluctuations and to increase the sensitivity of the receiver to detect a very weak radio signal transponded to the Earth 
from the spacecraft. On the other hand, technical limitations on the range of the transmitted frequency are not so 
restrictive and it can be changed smoothly in a very broad band according to a prescribed frequency modulation law. 
This law of modulation is chosen in such a way to ensure the receiving of the transponded signal from the spececraft 
exactly at the frequency v-i- It requires to know precisely the ephemerides of transmitter, observer, and spacecraft 
as well as the law of propagation of electromagnetic signal on its round-trip journey. Hence, one needs to know the 



Doppler shift 5v/v2 where 5v = vp, — V2- From equations (208), (209) we have 



5v vp dtn dti dto dl~2 

_ = _E_1= «_J:_^_^_l. (210) 

V2 V2 dip, dto dti dt2 
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As one can see from (210) there is no need to know explicitly the transformation between the proper time of the 
spacecarft, T%, and the barycentric time of the solar system, t\. This remark simplifies calculations. 

Accounting for relationship ( pOJ ) and expression (110) for the metric tensor yields at the point of emission 

,-1-1/2 

(211) 



dtp 
d% 



N 



N 



'l-vl)[l + 2Y ""V 1 ^ -4V-^ 



(1 - V • V a0 ) 2 



. .... , , \/l — v1 n r 0a — V a o • Too 

a— 1 / a— 1 V au 

where vo(to) is the barycentric velocity of emitter, v a o = v a (so) is the barycentric velocity of the a-th gravitating 
body, ro a = |ro a |, ^oa = xo(to) ~~ x a(so)j & n d So = to — ro a is the retarded time corresponding to the time of emission, 
to, of radio signal. 

Similar arguments give 
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a— 1 / a— 1 V aZ 

where v 2 (t 2 ) is the barycentric velocity of emitter, v a2 = v a (s 2 ) is the barycentric velocity of the a-th gravitating 
body, r 2a = |r 2a |, r 2a = x 2 (t 2 ) — x a (s 2 ), and s 2 = t 2 — r 2a is retarded time corresponding to the time, t 2 , of signal's 
reception. 

For up- and down- radio links the relationship ( pl| ) yields respectively 
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and 
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-Ba(s 2 ,si,t 2 ,k 2 ) 



Here the retarded time s\ comes out from the relation s± =ti — |xi — x (si)|, and 
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(215) 



|xi(ti)-xo(t )| ' ""' " 2 |x 2 (t 2 )-x 1 (t 1 )| ' 

are the unit vectors which define direction of propagation of transmitted and transponded radio signals respectively, 
and 



t x = to — ki • xq , 



to = ti — k 2 ■ xi 



(216) 
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2. Auxilary Partial Derivatives 



The relationships (p4|)-(p7|) allow to write down corresponding expressions for the retarded times sq, Si, and s 2 . 
One has to carefully dinstinguish between derivatives for the up- and down-radio links. For the transmitter-spacecraft 
up-radio link we have 
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These formulas must be used in equation (213). For the spacecraft-receiver down-radio link we obtain 
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These formulas must be used in equation (214). We point out that the meaning of the time derivative (217) is 



completely different from that of the time derivative (223) although they are calculated at one and the same point 
of transponding of the radio signal. At the first sight it may look surprising. However, if one remembers that the 



derivative ( 217 ) is calculated along the transmitter-spacecraft light path and that ( |223| ) along the spacecraft-receiver 
light path, which have opposite directions and different parameterizations, the difference becomes evident. 
The other set of time derivatives required in subsequent calculations reads as follows, 
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ata #21 ' 

where Rqi = |xo — xi| is the radial distance between emitter on the Earth and spacecraft, i? 2 i = |x 2 — Xi| is the radial 
distance between receiver on the Earth and spacecraft, and for the impact parameters hold ^ 1 = ki x (xi x ki), and 
t 2 = k 2 x (xi x k 2 ). 

Partial derivatives of functions B a (si, So,t*, k\) and B a (s 2l s±, t 2 , k\) can be found by making use of relationships 
©-(U). This yields 
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We have neglected in formulas (231),( |232| ) and (235), (236) all terms depending on accelerations of the solar system 
bodies. 
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3. Relativistic Effect for Doppler Measurement near Solar and Planetary Conjunctions 



The relationships ( 210 ) - (236) constitute the basic elements of the post-Minkowskian (Lorentz-covariant) Doppler 
tracking theory. They are sufficient to calculate the Doppler response for any conceivable relative configuartion of 
transmitter, spacecraft, and the solar system bodies. We shall consider in this section only the case when the spacecraft 
is beyond a massive solar system body like Sun, Jupiter or Saturn when the impact parameters of up- and down-radio 
links are small compared with distances from the body to transmitter, receiver, and spacecraft. We shall also rectrict 
ourselves to the consideration of gravitational shift of frequency only. Actually, this case is similar to gravitational 
lens. Thus, we neglect all terms of order m a /ro a , fn a /T-\_ a , rn a /r2a, fn a /Roi, vn a jRi\ as well as terms being quadratic 
with respect to the velocity v a . It is worthwhile to point out that the round trip travel time of the transmitted radio 
signal is much shorter than orbital period of any of the solar system body. For this reason, all functions with the 
retarded time argument entering the equations can be expanded around the time of transmission of the signal which 
is precisely determined by atomic clocks. Taking into account these remarks and making use of relationship ( |l69|) we 
obtain 

fSv\° s 2 n r i 4:Gm a 

I — = -K - ^ v o - pVi ■ a(£J , a (O = —3-2- £ Q , (237) 

\V2jgr C R R Cd a 

where v is velocity of transmitter, V! is velocity of spacecraft, v a is velocity of the a-th gravitating body deflecting 
trajectory of the emitted radio signal at the angle a 1 , d a = |£J is the length of the impact parameter of the light 
ray with respect to the a-th body, n is the distance between the transmitter and the light-deflecting body, ro is 



the distance between the spacecraft and light-deflecting body, and R = |xq — x — 1| ~ ro + r\. Formula (237) for 



Doppler shift by gravitational lensing depends on velocities of transmitter, spacecraft, and the ligh-deflecting body 



and generalizes that obtained independently by Bertotti & Giampieri [107] who considered only static gravitational 
lens. In case, of Doppler tracking observations of spacecraft in the field of Sun the difference between the two formulas 
is negligible, but the motion of the lensing body may be important in the case of Doppler observations of spacecrafts 
in the field of giant planets like Jupiter or Saturn. 

Approximate value of the Doppler shift is determined by the expression Sv/vz = 2a(v^/c) cosip, where a is the 
deflection angle of the light ray, v^ is velocity of the Earth, and tp is the angle between v^ and the impact parameter. 
For the Sun the deflection angle over the whole sky is not less than 1 milliarcsecond or ~ 4.85 • 10~ 9 radians. 
The relative velocity of the Earth with respect to speed of light is about 10~ 4 . These simple estimates applied to 
the Doppler shift's formula elucidate that the gravitational shift of frequency in Doppler tracking of interplanetary 
spacecraft caused by the Sun is not less than ~ 4.85 • 10~ 13 for any location of the spacecraft in the sky. If the path of 
the radio link grazes the Sun's surface the Doppler shift will be about 8.47- 10 -10 - a quantity which can be measured 
rather easily. The same kind of estimates gives for radio signals grazing the Jupiter and Saturn the Doppler shifts 
about 7.76 • 10~ 12 and 2.91 • 10~ 12 , respectively, which can be also measured in practice. 
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Our formalism for derivation of corresponding relationships for the description of high-precise Doppler tracking 
of interplanetary spacecrafts can be compared with approaches based on the post-Newtonian approximation scheme 



(see, e.g., [107], [127|). The advantage of the post-Minkowskian approach used in this paper is that it automatically 
accounts for all effects related to velocities of gravitating bodies through the expressions of the Lienard- Wiechert 
potentials. The post-Newtonian scheme makes calculations much longer and not so evident. 



4- Comparision of Two Mathematical Techniques for Calculation of the Doppler effect 

It is worthy from the methodological point of view to compare calculation of the Doppler effect in terms of frequency, 
used throughout the present paper, with that in terms of energy (see [p4| for definition) used, e.g. by Bertotti & 



Giampicri [107|. Let us introduce definitions of the 4- velocity of observer u a = u°(l,v' 1 ), the 4- velocity of source of 
light Uq = Uq(1,Vq), the 4- momenta of photon at the point of emission Kq = Kq(1, x z (to)) and the point of observation 
K, a = JC°(l,x l (t)), where u° = dt/dT, u° = dto/d%, JCq = dto/dXo and K.° = dt/dX with A and Ao being values of 
the affinc parameter along the light geodesic at the points of emission and observation. Then, using the definition of 



the Dopppler effect in terms of energy (275) it is not difficult to show that equation (275) can be recast into the form 

w /C°( 300 (i,x) + g0 i(t,x) [x*(t) +v i ] +ft J -(t,x)i i (t)t; i 
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fa 



u l^o\ 3oo(io,x ) +goi(tQ,x ) [i*(to)+^o] +9y( t o,xo)i'(*o)«'5 



(238) 



Calculation of time component K, of the 4- momentum of photon in ( 238) can be done if one knows the relationship of 
the affine parameter A along the light geodesic and coordinate time t. This is found by solution of the time component 
of the equation for the light geodesic (G = c = 1) 
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Using decomposition (H) of the metric tensor and parametrization (O) along the unperturbed light ray, equation 



(239) may be written 



ft_ 
dX 2 



-k a k^d t h al} - k a d T h 0a 
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where the constant vector k a — (1, fc l ) = (l,k), and the substitution for the unperturbed trajectory of light ray in 



h a is done after taking a partial derivative with respect to coordinate time t. Solution of equation ( 240 can be found 
by iterations using expansion 



A = E- 1 [t + F{t)] , 



(241) 



where E is the constant photon's energy at past null infinity measured by a fictitious observer being at rest, and the 
function T(t) is of the order 0(h a p). It is obtained by solution of the equation 
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d 2 T i 



dr 2 2 
Solving the differential equation (E42L) one finds 
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After examination of structure of integrands in the integrals of the expressions (244), (B45J) one notes that 
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A remarkable property of equality (p4q) is that the parameter t* is independent from the argument a of the integrand 



in (244), (245) and, for this reason, the derivative with respect to t* can be taken out of the sign of the integrals. It 
allows to transform, e.g., the integral (p44[) into the form 
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Using solution (|l(j) for h a p and relationship (|2^) relating total differentials of coordinate time a and retarded time £ 
one obtains 
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(248) 



t*+k-x a (C)-C 

where the upper limit s(r,t*) of integral in the right hand side is calculated by means of solution of equation, 
s + |k<7 + £ — x a (s)| = r + 1* , and depends on time r and instant of the closest approach t* considered as a parameter 



[ 148 1 . For the upper limit depends on t* the derivative d/dt* of the integral in square brackets is taken both from the 
integrand of the integral and its upper limit. It is possible to eliminate dependence of the upper limit of the integral 
on the parameter t* . It will be achieved if one takes time t as independent variable instead of r and finds the upper 



limit of the integral from the equation (83) as we have done previously while calculating the Doppler effect in terms 
of frequency. Such a procedure gives us 
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where the second term in the right hand side is a partial derivative of the upper limit of integral in (248) with respect 



to t* 1 149 1 and the upper limit s(t) of the last integral is treated as independent from t* |15C(|. Finally, one has 
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where the function C a (s) is displayed in (|7l|). Similar arguments give 
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Going back to the formula ( |238|) of the Doppler effect in terms of energy one can see that it can be factorized in three 
terms 
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1-k- v - v • S(t ) - k a h 0a (t ,x ) - k a v 3 h aj (t ,xo) 
Here H(r) is given in (J34J) and 3 (to) is obtained from ( p4| ) by means of calculation of all functions involved at the 
instant To- 

On the other hand, our previous result for calculation of the Doppler shift in terms of frequency obtained in section 
VI. C had the following form 



v _ dt 
v dt 



(257) 



Thus, in order to have an agreement with calcualtion of the Doppler shift in terms of energy one must prove that 

^-=S2-S 3 . (258) 

dt 



One can recast the product on the right hand side of J258|) accounting for equations ( {244] ) , (245), ( |250| ), ([251]) into 
the form 
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where the partial derivatives diB a p{r) and diB a p{T$) are calculated on the ground of equation (BQ). With equations 



( |259| ) - ( |263| ) it is straightforward to confirm the validity of equation (258) if one notes that up to the second order of 
the post-Minkowskian approximation scheme it holds 
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so that equation ( |25S| ) can be re-written as follows 
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It is easy to confirm that the numerator and denominator of ( |265[ ) coincide exactly with those of equation (pi]) used 
for calculation of the Doppler shift in terms of frequency and, for this reason, equation ( |258| ) is valid. This finalizes 
the proof of equivalence of using two different mathematichal techniques for calculation of the Doppler effect. 



In conclusion of this section we would like to point out that the method of calculation of integrals in formulas (244) 



fl245|) exposed in the sequence of equations (|246[) - (25C) significantly simplifies and reduces the amount of calculations 



which have been performed, e.g., in [151] for studying anisotropics of CMB radiation due to cosmic strings, where 
rather complicated transformations of variables were used for performing of the integrals under discussion. As we 
have shown in the present section such transformations are actually unnecessary. 

5. The Explicit Doppler Tracking Formula 

In view of practical applications it is useful to give the explicit formula for Doppler tracking of satellites. We shall 
derive it in the present section for one-way propagation of electromagnetic signals emitted from the point Xo at time 
£o and received at the point x at time t. The Doppler shift of the observed frequency v with respect to the emitted 



(>() 



frequency i/g is given by equation (252) which is to be transformed to separate the special relativistic Doppler effect 
from general relativistic corrections. Thus, we have 
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(266) 



where the first two factors out of four describe the special relativistic Doppler effect, and the next terms are general 
relativistic corrections. The unit vector k given at past null infinity relates to the unit vector K (see (|36j) for 
its definition) of the boundary value problem through the relationship ( p7\ ) which, for the particular case under 
discussion, reads 
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where R = |x — Xq|. 
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The explicit formulas for the functions a(r) and a(ro) are derived using (211) which leads to 
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We recall that Vo(io) is the barycentric velocity of emitter, v a o = v a (so) is the barycentric velocity of the a-th 
gravitating body at the instant sq, r$ a = |ro a |, i"o a = Xg(io) — x a(so)7 an d So = to — r^ a is the retarded time 
corresponding to the time of emission, to, of the radio signal. Besides this, v(i) is the barycentric velocity of receiver, 
v a = v a (s) is the barycentric velocity of the a-th gravitating body at the instant s, r a = |r a |, r a = x(i) — x a (s), and 
s = t — r a is the retarded time corresponding to the time of reception, t, of the radio signal. 



Omitting all terms in equation (74) for the integral C a depending on accelerations of the bodies' center-of-mass, 



and reducing similar terms, we obtain for the functions in the last factor of the basic relationship (266) the following 
explicit result 
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The formulas (266) -([271]) describe the Doppler shift of the radio signal transmitted from observer to spacecraft. 
The Doppler shift of the radio signal transponded back to the observer is described by a similar set of equations with 
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corresponding attachment of all quantities to the instant of the signal's reflection from the spacecraft and to the one 
of the signal's reception. In case of light grazing a gravitating body, the formula ( |266| ) gives, of course, the result 
shown already in ( |237| ). 

VIII. DISCUSSION 
A. Basic Results 

The long-standing problem of relativistic astrophysics and astrometry concerning propagation of electromagnetic 
signals in the weak but arbitrarily fast changing, time-dependent gravitational field of an astronomical N-body system 
is analytically solved in the present paper in the first post-Minkowskian approximation of General Relativity . The 
gravitational field, described by the perturbation h a p of the Minkowski metric tensor rj a p of the flat space-time, 
is presented in the form of the Lienard-Wiechert potentials and depends on coordinates, x a , (a = 1,2, ...,N), and 
velocities, v a , of the bodies taken at the retarded instants of time. There is no any restriction on the motion of 
the bodies except for that v a < c (speed of light). The relativistic equations of light propagation are integrated in 
the field of the Lienard- Wiechert potentials and their solution are found in algebraically closed form. Exact analytic 
expressions for the integrated time delay, the angle of light deflection, and the gravitational shift of electromagnetic 
frequency caused by the gravitational fields of arbitrary moving bodies are derived and all possible residual terms are 
shown explicitly. One can compare theoretical elegancy and completeness of the Lorentz-covariant formalism of the 
present paper with various approaches of other authors to the same problem of light propagation in time-dependent 



gravitational fields (sec, for example, [151] - [154]). 

The applications of the Lorentz-covariant theory of light propagation, developed in the present paper, to relativistic 
astrophysics and astrometry are as follows: 

• general theory of the Shapiro time delay in binary pulsars is developed and all corrections with respect to veloc- 
ities of pulsar and its companion to the standard logarithmic expression of the time delay in static gravitational 
field are found. Particular attention is paid to the terms being linear in velocities which generalize the formula 
for the Shapiro time delay which existed in the parameterized post-Keplerian formalism discussed by Damour & 
Taylor [173]. Lorentz-covariant post-Minkowskian approach to the time delay calculations is compared with the 
post-Newtonian approach the enigmatic efficiency of which remained puzzling for a long time, is fully explained 
both in terms of the analytic mathematical technique and in the visual language of Minkowski diagrams; 

• equation of gravitational lens, moving arbitrarily fast and possessing spin-dipole and quadrupole components, 
is derived. Gravitational shift of spectral lines of the lensed source of light is worked out and its influency on 
the anisotropy of cosmic microwave background radiation is discussed; 
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• the expression for the Shapiro time delay, caused by the solar system bodies, is re-analyzed to improve accuracy 
of pulsar timing data processing programs and of the consensus model of very long baseline interfcrometry; 

• relativistic deflection of light in the solar system gravitational field is obtained with accounting for all velocity- 
dependent terms in the first post-Minkowskian approximation. This result will be important in future space 
astrometric missions like GAIA (ESA), SIM (NASA), etc.; 

• theoretical formulation of the Doppler tracking of interplanetary spacecrafts is achieved at the level of residual 
terms of order 10 -16 . 

We could not elaborate in the present paper all possible aspects of the Lorentz-covariant approach to the problem of 
propagation of light rays in time-dependent gravitational fields of isolated astronomical systems. Some of the most 
important theoretical developments which can be done in future are outlined in the following section. 

B. Future Prospects 

The clear mathematical formulation of the Lorentz-covariant theory of light propagation in gravitational fields 
of arbitrary-moving bodies and the elegant method for solving related problems coming up in this framework and 
being based on the proper account for all retardation effects, delivers new fascinating opportunities for much deeper 
exploration of the following open problems of modern relativistic astrophysics and astrometry: 

• propagation of light rays in the field of arbitrary-moving bodies endowed with spin-dipole and quadrupole 
moments. This requires the knowledge of the expression for the singular tensor of energy-momentum of point- 
like particles with spin and quadrupole moments. Spin contribution to the tensor can be found, for example, in 
[H but the structure of the tensor with the quadrupole (and higher) multipole seems to be unknown. Solution 
of the given problem will admit a precise mathematical treatment of timing observations of a pulsar orbiting a 
Kerr black hole, as well as a unique interpretation of those X-ray and gamma-ray sources which are assumed to 
have a Kerr black hole at the center of their accretion disks related; 

• extension of the Lorentz-covariant theory presented in this paper to the event of strong gravitational fields. It 
will require finding solutions of the equations of light propagation in the second post-Minkowskian approximation 
of general relativity or another alternative theory of gravity. Here one can expect to find differences between 
predictions of two gravity theories which may be used for suggesting new observational tests of the theories. It is 
also interesting to note |3C]] that if the light-deflecting body or/and observer move too fast in a specific direction 
even the weak and hence linear gravitational field can become strong in a chosen coordinate frame. In such a 
case the linear post-Minkowskian approximation is not enough to give unambigious observational predictions of 
relativistic gravitational effects in propagation of light rays and a second iteration of the Einstein equations is 
required; 
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• elaboration of the formalism of the present paper on the case of polarized electromagnetic wave to calculate the 



rotation angle of the plane of polarization along the null geodesic path of the wave - the Skrotskii effect [155 



see also [151 1; 



• inclusion in the formalism of the given paper of relativistic effects of gravitational waves from localized sources 
like supernova explosion, massive binary black holes in nuclei of active galaxies, cataclysmic and ordinary binary 
stars in our galaxy, etc. The first decisive step towards the adequate interpretation of these gravitational wave 
effects has been done in our paper pj. However, a more involved technique is required to take into account 
motion of the sources of gravitational waves with respect to observer. We expect that new interesting effects 
may be found along this line; 

• calculation of response of space gravitational wave interferometers like LISA to the signals emitted by gravita- 
tionally induced oscillations of the Sun (so-called g- modes). The combined technique of this and our previous 
paper M is undoubtedly enough for getting the answer to that problem; 

• application of the formalism of the present paper to the case of small-angle scattering problem of fast-moving 
self-gravitating bodies and calculation of gravitational waveforms (cf., e.g., |p9[); 

• development of physically adequate, high-precision algorithms for data processing of observations of space as- 
trometric satellites and navigation systems like GPS as well as very long baseline interferometry. Practical 
necessity in such algorithms is strongly felt already today and will permanently grow following achievements in 
the rapid development of advanced space technologies. 

We could continue the list of subjects for future work. For example, we did not touch cosmological applications of 
the formalism of the present paper. This will require some modifications of equations of light propagation to account 
for cosmological expansion of the universe. No doubt, the interpretation of observations of anisotropy of cosmic 



microwave background radiation induced by, e.g. cosmic strings [151], can be made more theoretically adequate in 
the framework of the presented new scheme. 
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FIG. 1. Illustration of the light-ray's propagation history. The light ray is emitted at the instant of time to at the point 
xo and arrives at the point of observation x at the instant of time i. Light-deflecting bodies moves along accelerated world 
lines during the time of propagation of the light ray; their velocities at some intermediate instant of time are shown by black 
arrows. In the absence of the light-ray-deflecting bodies the light ray would propagate along an unperturbed path (dashed 
line) which is a straight line passing through the points of emission, xo, and observation, x. Direction of the unperturbed 
path is determined by the unit vector K = — (x — xo)/|x — xo|. In the presence of the light-ray-deflecting bodies the light 
ray propagates along the perturbed path (solid line). The perturbed trajectory of the light ray is bent and twisted due to 
the gravitoelectric (mass-induced) and gravitomagnetic (velocity-induced) fields of the bodies. The initial boundary condition 
for the equation of light propagation is determined by the unit vector k defined at past null infinity by means of a dynamical 
backward-in-time prolongation (dotted line) of the perturbed trajectory of light from the point of emission xo in such a way 
that the tangent vector of the prolongated trajectory coincides with that of the perturbed light-ray's trajectory at the point of 
emission. Relationship between unit vectors k and K includes relativistic bending of light and is given in the text by equation 

©■ 
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FIG. 2. Schematic illustration of the Shapiro time delay in a binary pulsar. The pulsar emits radio signal at the time to 
which reaches the observer at the time t. For calculation of the Shapiro time delay positions of the pulsar and its companion 
must be taken at the retarded instants of time so and s corresponding to those to and t. See also Fig. for further explanations. 
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FIG. 3. Schematic space-time diagram showing the relationship between positions of a photon taken at different instants of 
time to, ti,..., t§ (events 0,1,.. .,6 on the photon's null world line) and positions of the light-ray-deflecting bodies (marked by 
the black circles) taken at the instants of the retarded time corresponding to the instants to, t\,...,te. For simplicity only the 
two-body system is considered. The photon is deflected by the retarded gravitational field of the bodies expressed through the 
Lienard-Wiechert potentials. Also shown are positions of the bodies (marked by the unfilled circles) taken on the space-like 
hypersurfaces (dashed lines) of the time instants to, ti,..., te- As the photon approaches towards the system (events 0,1) it 
moves in the variable gravitational field of two bodies. After crossing the system (events 5,6) the gravitational field at the 
photon's position is "frozen" since the photon moves along the same light cone as the gravitational field propagates. The 
"freezing" of gravitational field takes place during propagation of the photon inside the system (events 2,3,4). Spatial positions 
of gravitating bodies taken at the retarded instants of time are very close to those taken at the hypersurfaces of constant time 
when photon moves near or inside the system. It explains why the post-Newtonian solution for the metric tensor can be applied 
in this situation as well as the post-Minkowskian one for calculation of the photon's propagation. Retarded and instanteneous 
spatial positions of the gravitating bodies are drastically different when photon is at large distance from the system (far outside 
the near zone). In this case only the post-Minkowskin retarded solution for the metric tensor can be applied for an adequate 
description of the gravitational perturbations of the photon's trajectory. 
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FIG. 4. Relative configuration of observer, source of light and a moving gravitational lens deflecting light rays which are 
emitted at the moment to at the point xo and received at the moment t at the point x. The lens moves along straight line 
with constant velocity from the retarded position x a (so) through that x a (s) and arrives to the point x a (i) at the moment of 
observation. Characteristic time of the process corresponds to the time of propagation of light from the point of emission up 
to the point of observation. 
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FIG. 5. The gravitational lens geometry for a moving lens M = J^ a=1 m a being at the distance r from the point of 
observation O with coordinates x l (t). A source of light S with coordinates xo(to) is at the distance R from O. Vector £ is the 
impact parameter of the unperturbed path of photon in the observer plane. Vector £ L denotes position of the center of mass 
of the lensing object in the lens plane. Vector r\ = BE is the observed image position of the background source of light S 
shifted in the lens plane from its true position by the gravitational field of the lens to the point E. Coordinates of the lens are 
X l (\*) = M~ x J2a=i m aXo(A*), and coordinates of the point E arc x l (A*) = x l (t) + s l (A* - t). 
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FIG. 6. Geodetic very long baseline interferometry measures delay r (the light travel time between points 2 and 3) in times 
of arrival of radio signal from a quasar at the first and second radio antennas, r = Ti — T\, located on the Earth' surface. 
Diurnal rotation and orbital motion of the Earth makes the delay to be dependent on time. This allows to determine the 
baseline b between two antennas, astrometric coordinates of the quasar, motion of the Earth's pole, parameters of precession 
and nutations, and many others. Modern data processing of VLBI observations is fully based on the relativistic conceptions 
and was supposed to be accurate up to 1 picosecond. 
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FIG. 7. Spacecraft Doppler tracking experiment in deep space. Radio signal is transmitted at the time to and at the point 
on the Earth along the unit vector ki . The radio signal reaches the spacecraft at the moment t\ and at the point 1 somewhere 
in the solar system and responds back to the Earth exactly at the same time ti along the unit vector k2 which has a different 
orientation from ki. The responded signal arrives at the reception point 2 on the Earth at the time £2- During the round-trip 
time of the radio signal the Earth rotates around its own axis and moves along the orbit. Hence, the barycentric position 
and velocity of the transmitter is different from the barycentric position and velocity of the receiver despite of that their 
topocentric positions on the Earth can coincide. When the impact parameter of the radio signal's trajectory is small the 
gravitational Doppler shift of the transmitted frequency with respect to the received frequency is estimated approximately as 
5v/v2 = 2q(vq/c) cos ip, where a is the deflection angle of the light ray, vq is velocity of the Earth, and <p is the angle between 
vq and the impact parameter. 
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